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Summary 

In  this  work  we  have  developed  (i)  a  preliminary  methodology  for  optimizing  multi-step  forging  process  and 
(ii)  the  corresponding  computational  tool  using  our  computing  substrate  PDESolve.  The  methodology  is 
based  on  the  solution  of  the  direct  and  sensitivity  problems  coupled  to  an  optimization  algorithm  and  is 
applicable  for  general  optimization-based  design  problems  in  mechanics.  In  particular,  it  can  be  used  for  the 
design  of  preform  shape  and  die  geometry  in  forging  processes.  Our  computing  substrate  PDESolve  is  an 
open  environment  programming  tool  based  on  the  object  oriented  language  C-I-+  and  has  specific  features  that 
allow  the  user  to  program  a  problem  at  a  high  level.  The  methodology  has  been  implemented  using  PDESolve 
and  will  be  applied  to  optimize  the  preform  shape  in  a  one-step  open  die  forging  operation  (upsetting).  The 
solution  of  this  example  as  well  as  others  in  mechanics  show  that  PDESolve  can  be  used  as  an  effective 
numerical  tool  for  solving  optimization-based  design  problems  in  engineering. 
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1  Introduction 

Forging  is  a  primary  deformation  process  that  plays  an  important  role  in  the  production  of  net-shape  or  near 
net-shape  components  for  aerospace  structural  applications.  During  metal  forging,  large  plastic  strains  are  im¬ 
parted  to  the  workpiece.  These  strains  cause  microstructural  changes  that  affect  the  material  response  during 
the  forging  process  itself  as  well  as  its  final  mechanical  properties.  High  strains  may  affect  the  workability  of 
the  workpiece  due  to  production  of  voids  and  fractures,  and  the  appearance  of  shear  bands  (flow  localization), 
all  of  which  will  result  in  a  defective  product.  On  the  other  hand,  material  properties  such  as  toughness  and 
strength,  as  well  as  the  formability  of  the  forged  part  during  secondary  processing  and  machining  are  very 
sensitive  to  the  developed  microstructure.  Moreover,  the  nonuniform  distribution  of  plastic  strains,  which  is 
typically  produced  by  most  deformation  processes,  may  induce  undesirable  residual  stresses  to  the  material 
which  will  affect  its  performance  in  service. 

Current  basic  research  at  the  Air  Force  is  emphasizing  the  development  of  fundamental  material  and 
process  models  to  address  many  of  these  aspects  so  as  to  improve  the  workability  of  the  material  as  well  as  the 
mechanical  properties  of  the  formed  part.  Forging  as  well  as  other  forming  processes  (rolling,  extrusion,  sheet 
forming)  are  being  investigated  for  a  number  of  materials  of  specific  interest  to  the  Air  Force,  such  as  various 
conventional  titanium  and  nickel-based  alloys.  Specific  aspects  of  the  research  focus  on  thermomechanical 
process  modeling,  microstructure  evolution  modeling  (grain  size,  grain  shape,  texture),  and  structure/property 
characterization  (constitutive  laws). 

Because  the  structural  performance  of  aerospace  systems  can  be  enhanced  by  a  judicious  tailoring  of  their 
shape  and  microstructure,  optimization  procedures  are  becoming  an  important  part  of  the  design  process.  In 
particular,  during  forging  the  final  geometry  and  state  (properties)  of  a  forged  part  depend  on  a  number  of 
processing  conditions,  such  as  forging  temperature,  initial  workpiece  geometry  (preform  shape),  die  geometry, 
lubrication  conditions,  loading  conditions  (ram  velocity),  as  well  as  the  type  of  material  being  forged.  De¬ 
signing  a  forging  process  to  meet  specific  requirements  of  performance  for  the  final  product  then  requires  the 
control/optimization  of  these  processing  variables, 

A  number  of  works  have  addressed  the  optimum  design  of  the  forging  process.  Kobayashi  et  al.  [21] 
introduce  the  so-called  backward  tracing  technique  for  preform  design.  Han  et  al  [17]  improved  on  this  method 
by  including  design  optimization  strategies  (sensitivity  analysis  and  specific  objective  functions)  and  solved  the 
design  of  intermediate  forging  die  shapes.  Both  of  these  methods  start  from  the  final  product  shape  and  trace 
back  the  loading  path  of  the  forging  process  to  determine  the  required  shape  of  the  preform  and  intermediate 
dies.  Recently,  sensitivity  analysis  and  optimization  techniques  have  been  coupled  with  the  solution  of  the 
direct  (forward)  problem  to  solve  the  preform  and  die  shape  design  problems  in  forging.  Applications  of  this 
approach  to  solve  preform  and  die  design  problems  are  given  by  Badrinarayanan  &  Zabaras  [2,  3]  for  one-step 
forging  and  by  Fourment  et  al  [10]  for  two-step  forging.  The  numerical  tool  used  in  all  these  studies  was  the 
finite  element  method. 

The  work  in  Phase  I,  which  will  be  continued  in  Phase  II,  is  aimed  at  developing  an  optimization-based 
forging  design  approach  to  determine  the  processing  variables  that  give  the  best  possible  compatibility  between 
final  product  shape  requirement  and  desired  material  state.  In  particular,  this  work  has  built  on  previous 
studies  [2,  3]  to  address  the  design  optimization  of  preform  shape  and  die  geometry.  Specifically,  in  this  Phase 
we  have  mainly  focussed  on 

•  The  development  of  an  optimization-based  design  methodology  for  multi-step  forging,  and 

•  The  implementation  of  the  methodology  using  our  computational  substrate  PDESolve  for  solving  2-D 
forging  design  problems. 

In  essence,  the  methodology  is  based  on  the  solution  of  the  direct  and  sensitivity  problems  defining  the  forging 
process  coupled  to  a  function  optimizer.  Because  the  present  version  of  PDESolve  is  based  on  the  finite 
difference  method,  the  governing  equations  for  the  direct  and  sensitivity  problems  have  been  implemented  in 
their  differential  (strong)  form.  Currently,  a  PDESolve  version  using  the  finite  element  method  (which  uses 
the  variational  or  weak  form  of  the  equations)  is  being  developed  and  will  be  used  as  the  numerical  tool  for 
Phase  II.  As  an  application  example  of  the  methodology,  the  PDESolve  code  has  been  used  to  solve  the 
preform  shape  design  for  open  die  forging  (upsetting),  modeled  as  a  2-D  problem.  This  application  uses  a 
hypoelastic- viscoplastic  model  for  the  material  description  and  assumes  sticking  friction  at  the  die-workpiece 
interface. 
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The  organization  of  the  report  is  as  follows.  Section  2  presents  a  brief  literature  review  on  sensitivity 
analysis  as  applied  to  solid  mechanics  problems.  Section  3  details  the  methodology  for  optimization-based 
design  of  multi-step  forging  together  with  the  formulation  of  the  equations  (in  strong  form)  for  solving  the 
direct  and  the  sensitivity  problems  in  nonlinear  solid  mechanics.  Section  4  presents  the  main  features  of  our 
computational  substrate  PDESolve  and  gives  some  application  problems  of  sensitivity  analysis  and  design 
optimization  in  solid  mechanics.  Here,  specific  examples  in  linear  elasticity  and  in  small  deformation  plasticity 
are  solved.  Section  5  gives  an  application  example  of  the  methodology  using  PDESolve  to  the  preform  design 
of  one-step  open  die  forging  (upsetting) .  Section  6  presents  a  road  map  for  an  extended  constitutive  framework 
to  be  included  in  Phase  II  and  Section  7  gives  the  conclusions  of  the  study. 

In  this  report,  we  use  direct  tensor  notation  [14]  to  develop  the  governing  equations.  For  this  purpose, 
vectors  and  second  order  tensors  axe  indicated  with  bold  face  lower  and  upper  case  letters,  respectively,  e.g. 
a  and  A;  and  fourth  order  tensors  are  denoted  by  calligraphic  bold  letters,  e.g.  C.  The  dyadic  product  of 
two  vectors  a  and  b  is  indicated  as  a  (8)  b.  Tensor  operations  between  two  second  order  tensors  A  and  B  are 
indicated  as  A-B  for  the  inner  product  (a  second  order  tensor),  A  (8)  B  for  the  dyadic  product  (a  fourth  order 
tensor),  and  A:B  for  the  scalar  product  (a  scalar).  Contraction  operations  over  two  indices  between  a  fourth 
C  and  second  A  order  tensors  are  represented  as  C:A  (a  second  order  tensor). 

2  Sensitivity  Analysis:  The  Basis  for  Optimization-Based  Design 
in  Solid  Mechanics 

Sensitivity  analysis  refers  to  methods  for  computing  the  rate  of  change  of  response  quantities,  such  as  displace¬ 
ments  and  stresses,  with  respect  to  changes  of  the  design  parameters.  This  first  derivative  information  (design 
sensitivity  coefficients)  is  used  in  gradient-based  design  optimization  algorithms  to  compute  the  gradient  of 
the  objective  function  and  constraints  with  respect  to  the  design  variables.  Typical  design  variables  in  solids 
mechanics  include  shape  (geometrical)  parameters,  sizing  parameters,  material  parameters,  and  boundary 
condition  parameters. 

Basically,  three  different  approaches  have  been  used  for  computing  design  sensitivity  coefficients  [18]:  The 
finite  difference  method  (FDM),  the  adjoint  structure  method  (ASM)  and  the  direct  differentiation  method 
(DDM).  In  the  FDM,  multiple  solutions  of  the  direct  problem  are  carried  out  with  different  values  of  a  single 
design  parameter  and  the  solutions  obtained  are  used  to  calculate  the  sensitivities.  The  number  of  evaluations 
required  depends  on  the  order  of  the  finite  difference  representation  used.  The  main  drawback  of  this  method 
is  the  critical  dependence  of  the  calculated  sensitivity  upon  the  step  size  of  the  design  variable.  Because  this 
method  is  costly  and  can  introduce  large  errors,  it  is  seldom  used  in  practice. 

The  ASM  is  based  on  formulating  an  auxiliary  (adjoint)  boundary  value  problem  for  the  sensitivity  co¬ 
efficients.  This  problem  is  constructed  by  imposing  certain  initial  stresses  or  strains  and  a  set  of  boundary 
conditions  defined  according  to  the  constraint  functional  (objective  function)  whose  design  sensitivity  is  de¬ 
sired.  The  solution  of  the  adjoint  structure  permits  then  the  explicit  evaluation  of  the  required  sensitivities 
for  the  direct  problem.  This  approach  is  very  effective  for  systems  defined  by  linear  field  equations,  such  as 
linear  elasticity  and  structural  mechanics  [18].  However,  it  is  not  suitable  for  history  dependent  problems 
since  the  adjoint  system  can  only  be  solved  after  the  analysis  of  the  direct  problem  has  been  completed  [30]. 
Computationally,  this  means  that  the  data  calculated  during  the  direct  analysis  need  to  be  either  saved  or 
recalculated  for  sensitivity  analysis.  This  is  quite  inefficient  for  computer  implementation.  In  the  DDM,  on 
the  other  hand,  the  sensitivity  equations  are  derived  by  taking  the  variations  of  the  continuum  field  equations 
of  the  direct  problem  (expressed  in  either  differential  or  variational  form)  with  respect  to  small  changes  of  the 
design  parameters.  The  resulting  equations  axe  linear  in  the  sensitivity  of  the  response  variables.  A  special 
feature  of  this  method  is  that,  for  history  dependent  problems,  the  sensitivity  equations  can  be  solved  once 
the  solution  of  the  direct  problem  for  a  particular  load  step  is  obtained  [30].  Hence,  the  direct  and  sensitivity 
problems  can  be  solved  in  parallel,  making  the  computer  implementation  of  the  method  quite  straightforward. 
For  this  reason,  the  DDM  is  recommended  for  highly  nonlinear  problems  with  a  path  dependence  response  [30], 
such  as  the  large  deformation  of  inelastic  solids. 

Sensitivity  analysis  and  optimization  for  linear  problems  in  mechanics  (linear  elasticity,  structural  me¬ 
chanics)  are  well  established  now.  A  number  of  excellent  books  [16,  18,  28]  have  been  published  that  gives  a 
rigorous  mathematical  treatment  of  sensitivity  analysis  for  linear  systems  [28]  as  well  as  its  application  to  solve 
engineering  problems  [16,  18].  In  this  area,  the  above  three  approaches  for  computing  sensitivities  together 
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with  optimization  techniques  have  been  widely  applied  for  engineering  design  using  the  finite  element  method 
and  the  boundary  element  method.  Particular  applications  are  in  the  automobile  and  aorospace  industries. 
The  reader  is  referred  to  a  recent  comprehensive  survey  paper  [8]  for  details. 

Sensitivity  analysis  and  optimization  techniques  as  applied  to  nonlinear  mechanics  have  recently  attracted 
the  attention  of  researches.  Quite  a  few  publications  deal  with  the  derivation  of  sensitivity  equations  for 
problems  involving  geometric  and  material  nonlinearities  in  structures  and  continua  [15,  20,  25,  30,  31,  32,  33]. 
Although  both  the  ASM  and  the  DDM  have  been  used  in  most  of  these  studies  [20,  30],  the  DDM  is  typically 
recommended  and  used  for  computer  implementation  with  the  finite  element  method  [1]  and  the  boundary 
element  method  [25].  Both  total  [30]  and  updated  Lagrangian  [20]  descriptions  have  been  presented.  Of 
particular  interest  for  the  present  work  is  the  study  reported  in  [2,  3]  which  deals  with  the  application  of 
sensitivity  angJysis  and  optimization  techniques  to  the  design  of  preform  shape  and  die  geometries  for  one- 
step  forging  processes.  This  study  will  be  used  as  a  guide  for  the  present  development. 

In  the  present  work,  we  mainly  focus  on  shape  design  sensitivity  analysis  using  the  DDM  for  nonlinear 
solid  mechanics  problems.  Our  particular  interest  is  in  modeling  the  finite  deformation  of  an  elcisto- viscoplastic 
material  as  a  means  to  analyze  metal  forging  processes.  In  this  work,  we  will  couple  the  solution  of  the  direct 
and  sensitivity  problems  with  an  optimizer  to  develop  an  optimization-based  design  approach  of  multi-step 
forging  processes. 

3  A  General  Framework  For  Optimization-Based  Design  Of  Form¬ 
ing  Processes:  Multi-Step  Forging 

3.1  General  Approach  for  Optimization-Based  Forging  Design 

The  methodology  for  solving  general  optimization-based  design  problems  requires  the  solution  of  the  governing 
equations  for  the  direct  and  sensitivity  problems,  coupled  to  a  function  optimization  algorithm.  For  the  specific 
case  of  the  forging  process,  the  developed  methodology  has  the  following  features: 

1.  The  direct  problem  is  formulated  as  the  large-deformation  analysis  of  an  elasto-viscoplastic  material. 
Under  isothermal  and  quasi-static  conditions,  this  deformation  process  is  governed  by  the  equilibrium 
equations,  the  kinematics  equations,  the  constitutive  relations  and  appropriate  initial  and  boundary 
conditions.  These  non-linear  differential  equations  are  written  using  an  updated  Lagrangian  description 
and  solved  numerically  using  an  incremental  strategy. 

2.  The  sensitivity  problem  is  formulated  using  the  DDM  [30].  In  this  method,  the  sensitivity  equations  are 
derived  by  considering  the  variations  of  the  governing  equations  of  the  direct  problem  with  respect  to 
small  changes  in  the  process  parameters  (design  variables).  As  mentioned  before,  this  method  has  proven 
adequate  for  highly-nonlinear  problems  with  a  path-dependent  material  response  [30].  The  sensitivity 
equations,  which  are  linear,  are  written  using  a  total  Lagrangian  description.  Their  numerical  solution 
is  also  carried  out  incrementally  and  in  parallel  to  the  solution  of  the  direct  problem. 

3.  The  minimization  of  the  objective  function  is  carried  out  using  standard  optimization  procedures  based 
on  either  a  gradient-based  method  or  a  quasi-Newton-based  method  [12]. 

In  this  work  (Phase  I),  we  focus  on  preform  (initial  workpiece)  and  die  shape  design  problems.  In  this 
context,  a  typical  shape  optimization  problem  can  be  stated  as: 

min  (t>{(3),  (3eR^  (1) 

subject  to  hj{f3)  =  0,  for  j  =  l,...,m  (2) 

^?A:(/3)>0,  for  k  =  l,...,p  (3) 

where  <l){f3)  is  the  objective  function,  /3  is  the  shape  design  variable  vector  with  n  components,  and  hj{f3)  and 
Qkif^)  are  equality  and  inequality  constraints,  respectively.  A  general  form  for  the  objective  and  constraints 
functions  in  non-linear  mechanics  can  be  represented  by  the  performance  functional 

m)=  f  /(u(/3),E(/3),T()3))dy 
Jv 


(4) 
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Figure  1:  A  two-step  forging  process. 


where  u,  E  and  T  represent  displacement,  strain  and  stress  tensor  fields,  respectively,  /  is  a  known  function 
of  its  arguments,  and  V  is  the  current  volume  of  the  body.  Here,  E  and  T  represent  some  general  strain  and 
stress  measures,  respectively. 

Standard  optimization  techniques  require  information  about  the  objective  function  and  its  derivatives.  In 
the  particular  case  of  shape  design  optimization,  the  derivatives  of  (j)  are  computed  by  taking  the  shape  design 
variations  of  (f)  with  respect  to  the  design  variables  /3.  Using  the  notation  <J(')  =  represent  the 

variation  (sensitivity)  of  the  quantity  (•)  with  respect  to  the  design  variation  can  be  written  as 


where 


S^=  {SfdV  +  fSdV) 

(5) 

Jv 

(6) 

SdV  =  tr{SL)dV  , 

(7) 

Here,  x  is  the  position  vector  of  a  particular  material  point  in  the  body,  and  SL  is  the  variation  of  the  design 
velocity  gradient  [3,  9].  The  design  velocity  represents  the  rate  of  change  of  the  position  vector  x  due  to 
changes  in  the  design  parameters,  i.e.,  v  =  ^.  Evaluating  the  functions  ^  and  5$  requires  the  solution  of 

the  governing  equations  for  the  direct  and  sensitivity  problems  to  compute  the  quantities  (u,  E,  T)  and  their 
corresponding  sensitivities  ((5u,  JE,(5T),  respectively. 

For  multi-step  forging,  the  process  design  is  formulated  as  the  design  of  the  preform  and  the  design  of  the 
dies  at  each  stage  of  the  forging  operation.  To  illustrate  the  procedure  to  be  followed  in  solving  the  multi-step 
forging  optimization  problem,  we  will  consider  the  two-step  forging  shown  in  fig.  1  modeled  as  a  2-D  problem. 
Here,  the  dies  are  assumed  to  be  rigid.  In  this  problem,  we  want  to  optimize  the  geometry  of  the  preform  and 
the  intermediate  dies  in  order  to  get  a  product  with  a  specified  geometry.  To  simplify  the  explanation,  we  will 
assume  that  the  geometry  of  the  initial  die  is  fixed.  The  particular  objective  function  for  this  optimization 
problem  can  be  written  as 

^  =  [  {n{f3,,M-n)^dv  (8) 

'  Jv 

where  H  is  the  computed  shape  from  process  simulation  and  is  an  implicit  function  of  the  response  field 
variables.  Cl  is  the  desired  shape,  and  (3^  and  13^  are  the  design  parameters  defining  the  preform  and  final  die 
geometry,  respectively.  Using  a  parametric  representation  for  these  shapes  such  as  Bezier  curves,  the  design 
parameters  can  be  identified  as  the  control  points  of  the  curves.  The  optimization  algorithm  will  then  require 
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F„ 


Figure  2:  A  general  step  of  a  forging  process. 


information  about  the  function  $  and  its  variation  (sensitivity)  with  respect  to  the  design  parameters  I3p  and 
/3d,  i.e., 

=  (9) 

To  compute  these  quantities,  we  need  to  solve  the  direct  and  sensitivity  problems  that  describe  the  large  strain 
inelastic  deformation  of  a  body.  The  body’s  deformation  will  be  described  using  the  updated  Lagrangian 
scheme  for  the  direct  problem  and  the  total  Lagrangian  scheme  for  the  sensitivity  problem. 

In  what  follows,  we  detail  the  framework  for  the  analysis  of  the  direct  and  sensitivity  problems  in  finite- 
deformation  elasto-viscoplasticity.  Then,  we  present  the  methodology  for  the  design  of  the  two-step  forging 
process  presented  in  fig.  1.  We  note  here  that  the  governing  equations  for  the  direct  and  sensitivity  problems 
will  be  presented  in  their  differential  form,  which  is  suitable  for  the  finite  difference  method  (Phase  I).  The 
variational  formulation  of  these  equations  to  be  used  with  the  finite  element  method  (Phase  II)  is  given  in  the 
appendix. 

3*2  Formulation  of  the  Direct  Problem 

The  direct  problem  consists  of  finding  the  deformation  history  and  material  state  of  a  body  subjected  to  the 
action  of  known  external  forces.  In  large  deformation  problems,  such  as  forging,  the  configuration  of  the  body 
changes  due  to  loading.  Figure  2  shows  a  general  step  of  a  forging  process.  Here,  we  denote  the  configurations 
of  the  body  in  the  undeformed  state  and  in  the  deformed  state  at  times  and  tiv+i  as  Bq^  Bn  and  Bjv+i, 
respectively.  We  assume  that  the  configurations  Bo  and  Bn  are  known  and  that  configuration  Bjv+i  needs 
to  be  determined.  A  particular  material  point  in  these  configurations  will  be  located  with  respect  to  a  fixed 
spatial  reference  frame  by  the  vector  positions  xo,Xiv,x/s^+i,  respectively. 

As  mentioned  above,  the  direct  problem  is  formulated  using  the  updated  Lagrangian  scheme.  Hence,  the 
equilibrium  equations  and  boundary  conditions  to  be  solved  for  determining  the  unknown  configuration  Bn+i 
are  written  relative  to  the  known  configuration  Bn  SiS 


Vtv  •  P  +  Pivb 

=  0 

in  Vn 

(10) 

P 

=  tN 

on 

(11) 

u 

=  u 

on  A^ 

(12) 

where  Vjv  =  9^57,  P  is  the  first  Piola-Kirchoff  stress  in  configuration  Byv+i  as  referred  to  the  configuration 
Bn  (P  =  Pjv^^)>  b  is  the  body  force  per  unit  mass,  pN  is  the  mass  density  in  configuration  Bn^  is  the  unit 
normal  to  the  surface  u  is  the  displacement  field,  and  tjv  and  u  are  prescribed  tractions  and  displacements 
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on  the  surfaces  and  respectively.  The  first  Piola-Kirchhoff  stress  is  expressed  as  (following  Gurtin’s 
convention  [14]) 

P  =  det{F)a  •  =  F  •  S  (13) 

where  F  is  the  relative  deformation  gradient  (F  =  Fj^"*"^),  a  is  the  Cauchy  stress  and  S  is  the  second  Piola- 
Kirchhoff  stress  in  configuration  Bat+i  as  referred  to  the  configuration  Bn  (S  =  We  note  here  that 

the  contact  between  the  die  and  the  workpiece  during  forging  is  introduced  in  eqs.(10)-(12)  as  a  traction-type 
boundary  condition. 

The  elasto-viscoplastic  constitutive  material  response  can  be  modeled  using  a  phenomenologically-based 
(continuum  plasticity  [23]  or  micromechanically-based  (polycrystal  plasticity  [5])  approach.  In  the  present 
work,  we  use  models  based  on  the  first  approach.  In  this  respect,  we  use  the  following  continuum  hypoelastic- 
viscoplastic  description  [22] 


a  =  C:(D-DP) 

(14) 

11 

Q 

(15) 

^  =  /(o-,s) 

s  =  g{a,  s) 

(16) 

(17) 

where 

?  =  <t_w<t  +  o-W 

(18) 

!l 

II 

(19) 

with  cr'  —  a  -  |tr(o')I.  In  these  equations,  D  and  are  the  rate  of  deformation  and  its  plastic  part, 

respectively,  W  is  the  continuum  spin,  ^  is  the  effective  plastic  strain  rate,  d  is  the  effective  stress,  s  is  the 
material  hardness  and  C  is  the  constant  elasticity  tensor  given  by  C  =  2pX  +  AI  (g)  I.  Here,  S  and  I  represent 
the  fourth-order  and  second-order  identity  tensors,  and  are  Lame’s  constants.  Also,  |1  er'  |1=  (or'io*')^/^. 

From  a  numerical  point  of  view,  the  time  integration  of  the  constitutive  equations  is  carried  out  using  an 
objective  integration  scheme  [34]  .  In  this  procedure,  the  constitutive  equations  (14)-(17)  are  transformed 
(pull  back)  using  the  rotation  tensor  R  obtained  from  F  =  R  •  U,  where  U  is  the  right  Cauchy-Green  strain 
tensor.  The  transformed  constitutive  equations  are  [22] 

^  =  C:(D-DP) 

(20) 

D-  = 

(21) 

^ 

s  =  g{a,  s) 

(22) 

(23) 

where  a  =  CT  R,  D  =  RT  .D  R,  and  similarly  for  and  N.  The  equations  for  the  rotated-neutralized 
tensors  a,  D,  and  N,  Eqs.(20)-(23),  are  then  integrated  using  an  implicit  integration  procedure  (backward 
Euler)  [22,  34]. 


The  numerical  solution  of  Eqs.(10)-(12)  and  Eqs.(20)-(23)  is  based  on  an  incremental,  displacement-based 
strategy.  For  this  purpose,  the  specified  time  interval  is  divided  into  a  certain  number  of  time  steps.  At 
each  time  step,  two  coupled  iterative  procedures  are  defined:  the  equilibrium  iterations  and  the  constitutive 
iterations. 

The  equilibrium  iterations  are  carried  out  to  compute  the  displacement  field  u.  Here,  we  use  the  differential 
fprm  of  these  equations,  written  in  a  suitable  form  for  algorithm  development.  This  form  is  obtained  as  follows. 
Using  the  updated  Lagrangian  nature  of  the  formulation,  we  can  write 

=  +  AF  =  I  +  AF 

=  +  AS  =  o-Af  +  AS 


F 

S 


(24) 

(25) 
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Hence,  P,  eq.(13),  can  be  expressed  as 

p  =  =  o-N  +  AS  -f  AF  *  ctn  +  AF  •  AS  (26) 

Neglecting  the  higher  order  term  AF  •  AS  in  Eq.(26),  we  can  write  Eqs(10)-(12)  as 

Vjv*(AS  + AF-^at)  =  <tn  inVN  (27) 

(AS  + AF-orjv) -utv  =  iN-cTN  -riN  onA%  (28) 

Au  =  Au  on  (29) 


Using  the  definition  of  the  second  Piola-Kirchhoff  stress  S  =  det(F)F”^  «  a  •  F”^,  the  pull  forward  of  the 
Cauchy  stress  cr  =  R  •  o'  •  R^,  the  polar  decomposition  of  the  relative  deformation  gradient  F  =  R  •  U,  and 
the  updated  Lagrangian  assumption,  one  can  show  that 


AS  =  Ad*  -  AU  •  cr  -  O'  •  AU  +  tr{AF)cr 
AU  =  (iF),  ,  AF=^ 


(30) 

(31) 


where  the  subscript  s  denotes  the  symmetric  part  of  (•).  In  Eq.(30),  the  stress  increment  Act  is  related  to  the 
strain  increment  AE  =  Odt  =  2(U  -  I)  •  (U  + 1)“^  using  the  elasto- viscoplastic  constitutive  matrix  C 
(consistent  tangent  operator  [27]) 

Act  =  £:A(AE)  (32) 

where  A(AE)  =  (AF),  =  AU  and 


C 

A 


2/iX  +  (/c  —  -A)I  <8i  I  —  2/i(77  —  c)N  ®  N 
o 


w  . 


b2 

CLib2  “}"  £12^1 


(33) 

(34) 


In  these  equations,  k  is  the  elastic  bulk  modulus,  a*  is  a  trial  effective  stress  (terminology  related  to  the 
constitutive  integration  scheme),  and  the  coefficients  01,02,61,62  depend  on  the  constitutive  functions  /  and 
g,  Eqs.(16)-(17), 

Using  Eqs.(30)-(34),  we  can  formulate  the  boundary  value  problem  (27)-(29)  in  terms  of  the  incremental 
displacement  Au  and  then  solve  the  resulting  equations  using  the  finite  difference  method.  The  current  version 
of  our  computational  substrate  PDESolve  will  be  used  for  this  purpose. 

The  constitutive  iterations,  on  the  other  hand,  are  carried  out  to  update  the  Cauchy  stress  cr  and  the 
material  hardness  s.  For  this  purpose,  the  bar-transformed  constitutive  equations  are  integrated  using  the 
backward  Euler  scheme  (radial  return  method) .  Details  of  the  algorithm  can  be  found  elsewhere  [22]  and  will 
not  be  presented  here.  In  the  context  of  software  development,  the  constitutive  integration  procedure  will 
be  added  to  PDESolve  as  a  particular  0+4*  class  and  then  coupled  to  the  solution  of  Eqs.(27)-(29)  as  an 
integral  part  of  the  solution  procedure. 


3.3  Formulation  of  the  Sensitivity  Problem 

The  formulation  for  the  sensitivity  problem  is  based  on  the  total  Lagrangian  scheme.  This  particular  scheme 
for  the  sensitivity  problem  avoids  difficulties  in  distinguishing  at  the  current  stage  of  the  process  between 
the  initial  shape  design  variation  and  those  due  to  the  deformation  process  itself.  Hence,  the  undeformed 
configuration  Bq  is  taken  as  the  reference  one  to  write  the  equilibrium  equations  and  the  boundary  conditions. 

In  the  present  work,  we  use  the  direct  differentiation  method  to  compute  the  sensitivity  of  the  equilibrium 
equations  and  the  boundary  conditions.  As  before,  the  notation  S{-)  will  be  used  to  represent  the  variation 
(sensitivity)  of  the  quantity  (•)  with  respect  to  the  design  variables  13^  and  (3^.  Taking  the  variation  of 
the  equilibrium  equations  and  boundary  conditions  written  in  the  configuration  Bn-^i  as  referred  to  the 
configuration  Bo,  we  obtain  the  following  sensitivity  equations 


Vo'(5P  ~  ^oP'AEq  +  —  0 

(<5P  -  P-<SLj)-no  =  rfto  +  [(rfF„-F-^):(I-n®n)- 

in  Vo 

(35) 

(F„-5Lo  •F"^):(n®n)]to 

on  Aq 

(36) 

*45 

II 

onAg 

(37) 
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where  Vo  =  P  is  the  first  Piola-Kirchhoff  stress  in  configuration  Bn+i  as  referred  to  configuration  Bo 
(P  =  P^^+i),  po  is  the  mass  density  in  configuration  no  is  the  unit  normal  to  the  surface  Aq,  to  is 
the  prescribed  traction  as  referred  to  configuration  Bq  (to  =  ^^t,  t  is  the  prescribed  traction  at  5),  u  is 
the  prescribed  displacement  and  is  the  total  deformation  gradient.  The  unit  normal  n  in  the  deformed 
configuration  Bn-\-i  is  related  to  no  by  n  =  det{Fu)^^^Z^  •  no.  Also,  in  these  equations,  SLq  represents  the 
variation  of  the  design  velocity  gradient  in  the  undeformed  configuration  Bq  and  is  given  by 

JLo  =  <5Ffl.F«'  =  ^  (38) 

where  Fr  =  is  the  mapping  of  material  points  from  a  reference  domain  Br  to  the  undeformed  configuration 
Bo-  This  reference  domain  shown  in  fig.  2  is  assumed  to  be  independent  of  both  the  design  parameters  and  the 
deformation  history  of  the  workpiece  (a  unit  square  domain,  for  example)  and  is  typically  used  when  dealing 
with  preform  shape  design  problems  [2].  In  eq.  (38),  xo  defines  the  initial  geometry  of  the  preform  (see  fig.  1) 
whose  shape  is  specified  by  the  design  parameters  /3p.  Hence,  ^Lo  is  a  function  of  these  design  variables.  If 
the  reference  domain  Br  coincides  with  the  undeformed  configuration  Bq,  then  Fr  =  1.  This  gives  5Lo  =  0. 
With  this  value,  the  sensitivity  boundary  value  problem,  eqs.(35)“(37),  reduces  to 


Vo*^P  "I"  Spoh  —  0 

in  Vo 

(39) 

JP  •  no  =  Jto  +  (JF„  •  F“^):(I  -  n  ®  n)to 

on 

(40) 

(Su  =  Su 

on  AS 

(41) 

Equations  (39)-(40)  are  mainly  used  when  dealing  with  die  design  problems  while  keeping  the  preform  shape 
fixed  [3]. 

Using  the  definition  of  P  =  det^Fu)^^  •  and  the  relations  cr  =  R •  o'  •  R^,  F  =  R-  U,  Ft^  =  F  -  Fn, 


one  can  show  that  [3] 

SP  =  detCFr,)[triSFn-FZ^)(T  +  6a-a-iSF^-F-Y]-'Pu'^  (42) 

6(t  =  R(JaR^  +  5RR^o--<T-5RR^  (43) 

JR  •  R’’  =  JF  •  F-^  -  R  -  JU  •  F-^  (44) 

JU  =  [U-^ -(F^-JF)*]*  (45) 

JF  =  6Fu-FZ^  -F^-F];,^  -SFn-FJ/-  (46) 


where  ^F^  is  related  to  the  displacement  variation  Sn  by 

SFu  =  ^-(F„-I)-(SLo  (47) 

This  formulation  is  complemented  with  the  sensitivity  equations  for  the  constitutive  model  which  relates 
the  stress  variation  Sa  to  the  strain  variation  5AE.  In  this  work,  these  equations  are  derived  from  the 
integrated  form  of  the  constitutive  equations  (20)-(23).  Omitting  details  of  the  derivation,  one  obtains 


S&  =  C:SAE-^A:SaN’^SalSi  (48) 

aiSsjsf)  (49) 

02 

where  Scr^  and  S$n  represent  the  variation  of  the  Cauchy  stress  and  the  hardness  at  time  respectively, 
and 

A  =  7?X+^-^I(g)I-(77-c)N(8)N  (50) 

SAE  =  4(U  +  I)“^  •(5U-(U-hI)-i  (52) 


As  mentioned  before,  the  sensitivity  problem,  Eqs.(35)“(52),  is  linear  in  Su  and,  hence,  its  solution  does 
not  require  an  iterative  procedure.  Also,  note  that  the  solution  of  the  sensitivity  equations  is  carried  out  in 
parallel  with  the  solution  of  the  direct  problem. 

When  solving  preform  and  die  design  problems,  the  design  parameters  and  (3^  are  introduced  in  the 
solution  procedure  of  either  eqs.(35)-(36)  or  eqs.(39)“(40)  in  two  ways: 
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1.  In  the  variation  of  the  design  velocity  gradient  (JLq  which  is  a  function  of  the  preform  geometry,  and 

2.  In  the  traction  boundary  conditions  along  the  workpiece-die  interface  which  depends  on  the  die  geometry. 

For  the  first  case,  we  note  that  the  mapping  F/j,  defining  the  reference  domain  Br,  accounts  for  the  geometric 
representation  of  the  preform,  and  hence,  is  expressed  in  terms  of  the  design  parameters  /3p.  As  an  example 
of  this  mapping,  consider  the  following  transfinite  mapping  [11]  of  the  preform  to  a  unit  square  (see  fig.l), 

a^oi  =  mG(m) 

X02  =  md  (54) 


where  G{7j2)  =  S/3pi^i(r72))  with  gi{r]2)  as  bcisis  functions,  defines  the  parametric  representation  of  the  free 
boundary  of  the  workpiece,  /3pi  is  the  i-th  design  parameter,  and  2d  is  the  height  of  the  preform.  Prom  this 
mapping,  we  can  derive 


Fr  = 


dxp 

dr) 


dG  1 


and  hence 


XT?  ^^^xa 


0  d 


9k  ^ 


Therefore,  the  expression  for  5Lo  =  (JF/j  •  is 

<JLo  = 


9k  -  ^9k) 

0  0 


Sl3pk 


The  sensitivity  equations  (35)-(37)  will  then  have  to  be  solved  for  each  I3pk  to  obtain  the  sensitivity  of  the 
field  variables  to 

For  the  second  case,  consider  the  traction  boundary  condition  along  the  interface  workpiece-die.  Assuming 
a  Coulomb  friction  model,  we  can  write  for  the  applied  traction  t, 

t  =  — pn±/x/ps  (55) 


where  p  in  the  normal  traction  (pressure)  to  the  contact  surface,  /x/  is  the  friction  coefficient,  and  n  and  s  are 
unit  vectors  normal  to  and  along  the  interface  respectively.  The  specific  sign  in  front  of  p/  is  chosen  according 
to  the  direction  of  motion  of  the  workpiece  with  respect  the  die.  Taking  the  design  variation  of  eq.(55)  with 
respect  to  the  design  parameters  and  using  the  relations:  n-n  =  1,  n-s  =  0,  one  obtains 


(5t  =  —[Sp  ±pf  p((Jn*s)]n  -f  [—p{Sn^s)  ±  PfSp]s  (56) 

Using  a  parametric  representation  for  the  die  surface  of  the  form  X2  =  h{xi,0^),  one  can  show  that  [2] 


Sns  = 


h^iSh^2  ““  h^2^h^i 


(57) 


where  h^i  =  ^  and  Sh^i  is  its  corresponding  design  variation.  Note  that  the  changes  in  the  die  surface  alters 
only  the  unit  normal  n.  By  substitution  of  eq,(57)  into  eqs.(39)-(41),  we  obtain  an  expression  for  the  design 
variation  of  the  traction  along  the  interface  as  a  function  of  the  design  parameters  13^.  In  this  case,  then, 
eqs.(39)-(41)  need  to  be  solved  for  each  parameter  f3d  to  determine  the  sensitivity  of  the  field  variables  to  the 
die  shape. 


3.4  Methodology  for  Optimization-Based  Design  of  Multi-Step  Forging  Pro¬ 
cesses 

The  ideas  behind  the  methodology  will  be  developed  by  considering  the  shape  design  optimization  of  the 
two-step  forging  process  represented  schematically  in  fig.  1.  Here  the  purpose  is  to  compute  the  field  variables 
(u,P)  and  their  sensitivities  (<Su,  <5P)  with  respect  to  the  design  parameters  (/3p,/3d)  at  the  end  of  the  two-step 
process.  This  information  is  then  used  by  the  optimizer  to  find  an  optimal  design  of  the  preform  and  die 
geometries. 
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This  problem  is  solved  using  eqs.(10)-(12)  for  the  direct  problem  together  with  eqs.(35)-(37)  and  eqs.(39)- 
(41)  for  the  sensitivity  analysis.  Computing  the  sensitivity  field  ((5u,JP)  with  respect  to  0p  during  the 
two-step  forging  can  be  considered  as  a  preform  design  process.  Similarly,  obtaining  the  sensitivities  (^u,  (5P) 
with  respect  to  (second  forging  step)  can  be  seen  as  a  die  design  problem.  Note  also  that,  during  the  second 
forging  step,  the  total  deformation  gradient  in  eqs.(35)-(37)  will  be  the  composition  ■  Fu  \  and  hence 
the  deformation  gradient  F^  ^  needs  to  be  saved  from  the  first  forging  step.  In  essence,  the  methodology  to 
solve  this  problem  is  as  follows: 

1.  Guess  the  initial  values  of  the  preform  design  parameters  /3p. 

2.  Do  i  =  1,  #  stepsl  (first  forging  step) 

(a)  Solve  the  direct  problem  for  (u,P).  Use  eqs.(10)-(12). 

(b)  Solve  Tip  sensitivity  problems  for  the  variations  ((5u,  JP)/?^,  where  rip  is  the  number  of  preform 
design  parameters.  Use  eqs.(35)-(37). 

3.  At  end  of  first  forging  step,  the  sensitivities  (iu,  (SP)^p  will  be  used  as  initial  conditions  for  solving  the 
second  forging  step. 

4-  Guess  the  initial  values  of  the  die  design  parameters  or  use  previous  computed  value. 

5.  Do  i  =  1,  #  steps2  (second  forging  step) 

(a)  Solve  the  direct  problem  for  (u,P).  Use  eqs.(10)-(12). 

(b)  Solve  Tip  sensitivity  problems  for  the  variations  (Ju,  (5P)^p.  Use  eqs.(35)-(37). 

(c)  Solve  rid  sensitivity  problems  for  the  variations  (Ju,  (SP)/?^,  where  rid  is  the  number  of  die  design 
parameters.  Use  eqs.(39)-(41). 

6.  Use  an  optimization  algorithm  to  compute  new  values  of  (/3p,/3^).  Here,  the  objective  function  ^  is 
minimized  using  gradient  information 

7.  Check  convergence  of  the  objective  function: 

If  II  $  II <  Tolerance  =>  exit 

K  11  ^  lj>  Tolerance  go  to  step  2  with  the  new  values  of  13 p  and 

4  A  Computational  Tool  for  Solving  Optimization-Based  Design 
Problems  in  Mechanics:  PDESolve 

4.1  General  Features  of  PDESolve 

PDESolve  is  an  object-oriented  class  library  constructed  using  the  C++  language  and  designed  for  solving 
complex  partial  differential  equations  (PDE)  with  general  tensor  value  functions.  It  forms  the  basis  for  an 
open-standard  computer  language  for  scientific,  engineering,  and  financial  PDE  applications.  When  using 
PDESolve  the  user  specifies  the  problem  to  be  solved  in  a  high-level  form,  together  with  the  discretization 
methods  and  the  solution  procedures.  This  mathematical  level  of  abstraction  gives  the  user  a  fully  flexible 
programming  language  for  constructing  PDE  simulations  rapidly.  Besides,  the  object-oriented  architecture 
of  PDESolve  allows  a  user  to  incorporate  new  capabilities  into  an  existing  PDESolve  code  in  a  clean  and 
efficient  manner. 

The  basic  C++  classes  in  PDESolve  are: 

•  Dif  f  Op:  represents  differential  operators 

•  Expr:  builds  mathematical  expressions  including  PDEs 

•  SpatialDomain:  specifies  the  problem  domain  (geometry) 

•  BC:  prescribes  boundary  conditions 


BEAM  Technologies,  Inc,  PROPRIETARY 


12 


•  FDMesh:  specifies  the  finite  difference  mesh  (spatial  discretization) 

•  DiscExpr:  constructs  the  discrete  form  of  the  PDE  using  the  mesh  and  the  boundary  conditions 

•  Function:  represents  both  known  and  unknown  functions 

Differential  Operators  and  Expressions:  In  PDESolve,  the  partial  differential  operator  is  rep¬ 

resented  as: 

DiffOp  dxx(2,0,4); 

In  the  argument  list,  the  first  integer  specifies  the  order  of  the  differential  operator,  the  second  integer,  the 
index  of  the  independent  variable  (i.e.  x  is  taken  to  be  independent  variable  #  0),  and  the  third  integer, 
the  order  of  accuracy  desired  in  the  numerical  representation  of  the  differential  operator.  Currently,  up  to 
fourth-order  accuracy  is  supported  for  the  differential  operators.  Optionally,  a  user  defined  stencil  family  file 
may  be  provided.  The  syntax  for  defining  the  3-D  Laplace’s  equation  to  third-order  accuracy  is 

DiffOp  dxx(2,0,3); 

DiffOp  dyy(2,l,3) ; 

DiffOp  dzz(2,2,3); 

Expr  laplace  =  (dxx  +  dyy  +  dzz)  ♦  u; 

Geometry  Specification:  PDESolve  uses  a  multi-block  body-fitted  coordinate  representation  of  the  problem 
geometry.  This  representation  is  based  on  two  binary  operators,  (Cartesian  product)  and  “-1-”  (union), 
together  with  a  base  operator  for  creating  line  segments  SpatialDomain(<high>,<low>).  With  this  repre¬ 
sentation,  the  geometry  and  discretization  are  specified  unambiguously.  The  PDESolve  syntax  for  defining 
a  ID  domain  (0,1)  is 

SpatialDomain  geom(0.0,1.0) ; 

A  unit  cube  is  specified  as 

SpatialDomain  geom  =  SpatialDomain (0.0, 1.0) 

*  SpatialDomain(0.0, 1 .0) 

*  SpatialDomain(0 .0,1.0); 

More  complicated  shapes  can  be  generated  by  adding  together  components  such  as  in  the  following  definition 
of  an  “L”  shaped  geometric  domain: 

SpatialDomain  geom  =  SpatialDomain (0.0, 1.0)  *  SpatialDomain(0.0, 1 .0)  + 

SpatialDomain (1.0, 2.0)  *  SpatialDomain (0.0, 1.0)  + 

SpatialDomain (0.0, 1.0)  *  Spat ialDomain(l .0,2.0) ; 

PDESolve  can  be  used  to  solve  problems  on  body-fitted  meshes  around  complex  shapes  as  well.  This  is  done 
through  a  user-specified  transformation  map  that  could  be  analytical  or  numerical. 

Boundary  Conditions:  For  the  line  segment,  to  specify  u  =  0  at  the  left  boundary  and  to  leave  the 
right-boundary  value  unconstrained,  the  syntax  is: 

BC  boundaryCond(u“0,  Unc()); 

For  the  unit  cube,  boundary  conditions  can  be  of  form 


BC  be  =  BC(u==0,  u==l) 

*  BC(u==0,  Dy*u==0) 

*  BC(u”0,  Dz*u==0)  ; 

where  the  Dif  f  Op’s  Dy  and  Dz  are  used  to  specify  Neumann  boundary  conditions.  The  PDESolve  boundary 
condition  language  is  based  on  structural  parallelism  with  the  geometry  specification.  This  parallel  structure 
provides  a  well-defined  specification  of  the  boundary  conditions. 

DiscExpr:  Discretized  Version  of  PDE  Problem:  Once  the  PDE,  mesh  and  boundary  conditions  have 
been  specified,  the  discretized  problem  is  obtained  using  the  statement 
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DiscExpr  discprob(laplace,  grid,  u,  bcs) ; 

The  above  command  causes  the  discrete  equations  to  be  built  and  stored  as  a  DiscExpr  object  called  discprob. 
This  system  can  be  solved  for  u  at  the  mesh  points  simply  by  specifying 

discprob. solve 0 ; 

Simple  Example  in  PDESolve:  We  now  consider  a  simple  nonlinear  PDE  problem  to  illustrate  how  the 
above  C+-1-  classes  are  put  together  to  solve  a  PDE  problem.  We  solve 


d^u  du  ... 


(58) 


on  the  interval  [0, 1],  where  f{x)  =  x{x  -  10)  “  4,  with  the  boundary  conditions  being  u\o  ==  0,u|i  =  1.  The 
executable  PDESolve  code  for  this  problem  is: 


#include  "pdesolve.h” 

Real  fFunc (const  Coords  &X)  { 
return  X[0]*(X[0]  -  10.0)  -4}; 
void  mainO  { 

Spat ialDomain  domain (0,  1); 

FDMesh  grid(domain,FDGrid(ll)) ; 

Function  u; 

Function  f (1 , scalar ,fFunc) ; 

BC  bc(u==0,  u==l) ; 

DiffOp  dx(l,0,2); 

DiffOp  dxx(2,0,2); 

Expr  eq  =  dxx*u  +  u*dx*u  -  f ; 

DiscExpr  discExpr(eq,  grid,  be,  NewtonSolver(40,le-5)) ; 
u  =  discExpr. solve () ;  //Solve  for  u 
cout  «  u  «  endl; 

} 


Note  that  here  we  specified  that  the  non-linear  problem  be  solved  with  Newton’s  method.  The  linearization 
is  performed  automatically  (although  the  user  can  input  the  linearization  and  loop  “by  hand”)  with  the 
maximum  number  of  iterations  being  40  and  the  tolerance  being  10“^.  Any  other  nonlinear  solver  can  be 
used  by  suitably  specifying  it  in  the  last  argument  of  the  DiscExpr.  This  argument  is  optional  (with  Newton 
linearization  being  the  default)  and  is  not  present  for  linear  problems. 

If  we  instead  need  to  solve  the  PDE 

+  u—  =  f{x),  (59) 


dx^  dx 


the  only  line  that  needs  to  be  changed  is 


Expr  eq  =  dxxx*u  +  u*dx*u; 

remembering  to  include  the  definition  for  dxxx  and  to  provide  an  additional  boundary  condition: 

BC  bc(u”0  &&  dx*u==l,  u"l)  ; . 

This  will  specify  that  on  the  left  side  of  the  domain  u=0  and  dx*u=l. 

In  the  next  sections,  we  use  PDESolve  to  solve  optimization-based  design  problems  in  linear  elasticity  and 
small  deformation  plasticity.  We  will  apply  the  general  methodology  developed  for  optimization-based  design 
problems  to  both  cases  and  show  how  PDESolve  can  be  used  as  an  effective  numerical  tool  to  solve  this  kind 
of  problems.  As  mentioned  before,  we  will  use  the  differential  (strong)  form  of  the  governing  equations  which 
are  amenable  to  the  finite  difference  method  (current  version  of  PDESolve). 
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4.2  Design  Optimization  in  Linear  Elasticity 

4.2.1  Problem  Formulation 

For  linear  elastic  problems,  the  different  configurations  of  the  body  shown  in  fig.  2  will  coincide,  and  hence, 
Bo  =  Bn  =  Bn-^i  (xo  =  Xiv  -  xyv+i  =  x).  For  this  case,  the  strain  and  stress  measures  are  given  by  the  small 
elastic  strain  e  and  the  Cauchy  stress  cr.  Hence,  the  general  performance  functional,  Eq.(4),  can  be  written  as 

$(/3)=  f  /(u(/3),c(/3),tT(/3))dy  (60) 

Jv 

As  before,  for  optimizing  this  functional  we  need  to  solve  both  the  direct  problem  for  the  quantities  (u,  6,(r) 
and  the  sensitivity  problem  for  the  variations  (<5u,<5€,(5o*). 

For  an  isothermal,  quasi-static  deformation  process,  the  direct  problem  of  linear  elasticity  is  governed  by 
the  equilibrium  equation  and  boundary  conditions  written  as 


V*(T  -h  f  =  0  in  V 

(61) 

cr-n  =  t  on  Ao^ 

(62) 

u  =  u  on  Au 

(63) 

together  with  the  constitutive  material  response  given  by 

o-  =  C:€  ,  e  =  i(Vu+(Vuf) 

(64) 

In  these  equations,  f  is  the  body  force  per  unit  volume,  n  is  a  unit  normal,  t  and  u  are  prescribed  tractions 
and  displacements  on  the  surfaces  A^  and  Au,  respectively,  and  C  is  the  constant  elasticity  matrix.  Note  that 
these  equations  are  linear  and,  hence,  their  solution  can  be  obtained  in  one  step. 

The  sensitivity  problem  is  formulated  using  the  direct  differentiation  method.  Here  the  sensitivity  equations 
are  obtained  by  taking  the  variations  of  the  equations  governing  the  direct  problem,  eqs.(61)-(64),  to  give,  for 
the  equilibrium  equations  and  boundary  conditions 


V-(5(t  -  Vff;(5L'^  +  (5f  = 

0  in  y 

(65) 

{Sc  —  <r*<5L^)-n  +  (JLin  ®  n)t  = 

5t  on  Aa^ 

(66) 

Su  = 

Su  on  Au 

(67) 

and  for  the  elastic  constitutive  response 

Sa  =  CiSe  ,  Se  =  i[V<5u  +  (VJu)^  - 

A 

(Vu-<5L  +  (Vu-<JL)^)] 

(68) 

where  SL  =  SFr-Fj^^  =  In  this  case,  the  design  parameters  (3  are  introduced  in  the  problem  through  the 
mapping  Fr  (we  use  transfinite  transformations  for  this  purpose). 


4.2.2  Implementation  of  2-D  Linear  Elasticity  Problems  in  PDESolve 

The  foregoing  equations,  eqs.(61)-(68),  have  been  simplified  for  the  case  of  two-dimensional  (plane  stress  and 
plane  strain)  problems.  The  expanded  form  of  these  equations  is  given  by: 

Direct  Problem 
Equilibrium  equations, 

^11,1  +  0^12,2  + /l  =  0 

<^21,1  ^22,2  + /2  ==  0 

Boundary  conditions, 

<^11^1  +  cri2^2  = 

+  cr22Tl2  =  ^2 
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Constitutive  law, 

an  =  (2// +  A)eii -f  A€22 

^22  =  Aeii  -f  {2fji  +  A)c22 

ai2  =  2//ei2 

where  A  =  A  for  plane  strain,  and  A  =  for  plane  stress, 

Kinematics  relations, 

€ll  =  ui,i 
€22  =  ^^2,2 

1/  N 

€12  =  +  W2,l) 

Sensitivity  Problem 
Equilibrium  equations, 

San, i  ^^12,2  ~  ^  ai2^iSLi2  —  ai2^2^L22  ^fi  ”  ^ 

Sa21,l  H"  ^^“22, 2  “  <^2l,lSLn  ”  ^21,2SL21  ^  <y22,lSLi2  —  Cr22,2SL22  +  Sf2  =  0 

Boundary  conditions, 

Sail'll  Sai2Ti2  H“  (^^12^-1  —  an^2)[{SLn'^2  ”  6Li2Tii)7ii  -h  (<5i/2i^2  ^-^'22^1)^12]  “  ^^1 

(5(721^1  +  5cr22n2  +  (o‘22ni  —  a2lTl2)[{SLn'^2  SLi2Tli)ni  +  {SL21TI2  —  <^^22^11)^2]  —  St2 

Sui  =  Sui 

Su2  =  Su2 

Constitutive  law, 

San  —  (2/i  +  A)(Seii  “h  A5e22 
Sa22  —  AJcii  4”  (2/x  +  A)(5€22 
Sai2  =  2fiSei2 

Kinematics  relations, 

^6ii  =  —  ui^iSLn  ”  ‘^i,2^'^2i 

^€22  =  <^^2,2  ^2, 1^1/12  “  U2,2SL22 

S€i2  =  2{Sui^2  +  <^W2,l)  —  —iUi^iSLi2  +  Ui^2SL22  +  U2,lSLn  4*  '^2,2SL2\) 

The  implementation  of  these  equations  has  been  carried  out  in  PDESolve  using  existing  C+H-  classes. 
The  corresponding  modules  are  shown  below.  The  specified  boundary  conditions  presented  in  this  program 
correspond  to  the  case  of  a  thin  circular  plate  with  a  hole  (plane  stress)  loaded  along  the  outer  radius  (see 
Fig.  3).  For  this  problem,  SLij  =  is  computed  using  the  transfinite  bilinear  mapping, 

Xi  =  [a  +  7j2(6-o)]cos(^) 

X2  =  [a  +  r,2(6-a)]sin(^) 

where  a(=  /3i)  and  b{=  P2)  are  the  inner  radius  and  outer  radius  of  the  plate,  respectively.  This  program  has 
been  used  to  predict  the  sensitivity  of  the  radial  and  tangential  stresses,  Sarr  and  Saeo^  to  variations  of  the 
inner  radius  (/3i).  Comparison  of  computed  results  with  the  analytic  solution  is  shown  in  Fig.  4. 
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Circular  Plate  With  Hole 


<5rr 


Figure  3:  Circular  plate  with  a  hole.  Figure  4:  Stress  sensitivities  for  circular  plate  with  a  hole. 

void  Elasticity_Direct() 

{  U  :  displacement  vector,  (U0,U1) 

E  :  elastic  strain  tensor,  (E00,E11 ,E01) 

T  :  Cauchy  stress  tensor,  (T00,T11,T01) 

F  :  body  force  vector,  (F0,F1) 

Po  :  load  applied  to  outer  radius,  (PoO,Pol) 

M  :  ‘3x3  elasticity  matrix  (plane  stress)  *\ 

SpatialDomain  computDomain=Spat ialDomain (0 . , 1 . ) *Spat ialDomain (0 . , 1 . ) ; 

Function  fwdMap(2,Tensor (2, 1) ,fwdMapFunc) ,Function(2,Tensor(2,2) ,DfwdMapFunc) ) ; 

Function  invMap (2, Tensor (2,1) .invMapFunc) ; 

SpatialDomain  physDomain=SpatialDomain (computDomain ,Map  (f wdMap ,  invMap)  )  ; 

FDMesh  grid=FDMesh (physDomain , FDMeshSpec (nX) *FDMeshSpec (nY) ) ; 

DiffOp  Dx(l,0,4); 

DiffOp  Dy(l,l,4); 

Function  U(2,Tensor(2,l)) ; 

Function  E(2,Tensor(3,l)) ; 

Function  T(2,Tensor (3, 1) ) ; 

Function  F(2,Tensor (2,1) ,CFuncF) ; 

Function  Po (2, Tensor (2,1) ,CFuncPo) ; 

Function  M (2, Tensor (2,3) ,CFuncM) ; 

Expr  diff0perDirect=List(List(Dx,0.0,Dy) ,List (0,0,Dy ,Dx)) ; 

Expr  E=tremspose  (dif  f  OperDirect)  ♦!! ; 

Expr  T=M*E; 

Expr  equilEqns=diff OperDirect *T+F ; 

BC  bc=BC (T[0]*Normal(0)+T [2] ♦Normal (1) ==0.0  &&  U[1]==0.0, 

U[0]==0.0  &&  T[2]^Normal(0)+T[l]+Normal(l)==0.0) 

♦  BC(T[0] ♦Normal (0)+T[2]+Normal(l)==0.0  && 

T  [2]  ♦Normal  (0)  +T  [1]  ♦Normal  ( 1 )  ==0 , 0 , 

T  [0]  ♦Normal  (0)  +T  [2]  ♦Normal  (1 )  ==Po  [0]  && 

T  [2]  ♦Normal  (0)  +T  [1]  ♦Normal  (1)  ==Po  [1]  )  ; 

DiscExpr  elastDirect (equilEqns ,  grid,  U,  be); 

U  =  elastDirect. solve  0 ; 
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void  Elasticity_Sensitivity() 

{  \* 

vU  :  variation  of  displacement,  (vU0,vUl) 

vE  :  variation  of  strain,  (vE00,vEll,vE01) 

vT  ;  variation  of  stress,  (vT00,vTll,vT01) 

vF  :  variation  of  body  force,  (vF0,vFl) 

vPo  ;  vairiation  of  applied  traction,  (vPoO,vPol) 

vL  :  variation  of  the  design  velocity  gradient,  ((vL00,vL01) , (vL10,vLll)) 

M  :  3x3  elasticity  matrix  (plane  stress) 

*\ 

SpatialDomain  computDomain=SpatialDomain(0,l)*SpatialDomain(0, 1) ; 

Function  fwdMap (2, Tensor (2,1) ,fwdMapFunc,Function(2,Tensor(2,2) .DfwdMapFunc)) ; 
Function  invMap(2,Tensor(2,l) ,invMapFunc) ; 

SpatialDomain  physDomain (comput Domain ,  Map (fwdMap , invMap) ) ; 

FDMesh  grid=FDMesh (physDomain, FDMeshSpec (nX) *FDMeshSpec (nY) ) ; 

DiffOp  Dx(l,0,4); 

DiffOp  Dy(l,l,4) ; 

Function  vU(2,Tensor(2,l) ) ; 

Function  vE(2,Tensor (3, 1)) ; 

Function  vT(2,Tensor (3, 1) ) ; 

Function  vF (2, Tensor (2,1) ,CVFuncF) ; 

Function  vPo (2, Tensor (2,1) ,CVFuncPo) ; 

Function  vL (2, Tensor (2,2) ,CVFuncL) ; 

Function  M (2, Tensor (2, 3) ,CFuncM) ; 

Expr  operLl=List(List(vL[0]  [0]  ,0.0,vL[0]  [1])  ,List(0.0,vL[0]  [1]  ,vL[0]  [0] ))  ; 
Expr  operL2=List(List(vL[l]  [0]  ,0.0,vL[l]  [1] ), List (0.0, vL[l3  [l],vL[l]  [0])); 
Expr  dif f 0perDirect=List (List (Dx ,0.0, Dy) , List (0 . 0 , Dy , Dx) ) ; 

Expr  dif f 0perSensit=operLl*Dx+operL2*Dy ; 

Expr  vE=transpose (dif f OperDirect) ♦vU-tremspose (dif f OperSensit) *11 ; 

Expr  vT=M*vE; 

Expr  vEquilEqns=diff OperDirect ♦vT-diff OperSensit *T+vF ; 

BC  vBC=BC (vT [0] ♦Normal (0) +vT [2] ♦Normal (1) + (T [2] ♦Normal (0) -T [0] ♦Normal (1) ) ♦ 
((vL[0] [0] ♦Normal (l)-vL[0] [1] ♦Normal (0)) ♦Normal (0)+ 

(vL[l]  [0] ♦Normal (l)“vL[l] [1] ♦Normal (0)) ♦Normal (1)) ==0.0  && 
vU[l]==0.0, 
vU[0]==0.0  && 

vT  [2]  ♦Normal  (0)  +vT  [1]  ♦Normal  ( 1)  +  (T  [1]  ♦Normal  (0)  -T  [2]  ♦Normal  ( 1)  )  ♦ 

(  (vL  [0]  [0]  ♦Normal  (1)  -vL  [0]  [1]  ♦Normal  (0)  )  ♦Normal  (0)  + 

(vL[l]  [0] ♦Normal (l)-vL[l]  [1] ♦Normal (0)) ♦Normal (1))  ==0.0)) 

♦  BC  (vT  [0]  ♦Normal  (0)  +vT  [2]  ♦Normal  ( 1)  +  (T  [2]  ♦Normal  (0)  -T  [0]  ♦Normal  ( 1)  )  ♦ 

(  (vL  [0]  [0]  ♦Normal  (1)  -vL  [0]  [1]  ♦Normal  (0)  )  ♦Normal  (0)  + 

(vL  [1]  [0]  ♦Normal  (1)  -vL  [1]  [1]  ♦Normal  (0)  )  ♦Normal  (1) )  ==0 . 0  && 
vT  [2]  ♦Normal  (0)  +vT  [1]  ♦Normal  ( 1)  +  (T  [1]  ♦Normal  (0)  -T  [2]  ♦Normal  ( 1)  )  ♦ 
((vL[0]  [0]  ♦Normal (l)“vL[0]  [1]  ♦Normal  (0))  ♦Normal  (0)  + 

(vL[l]  [0] ♦Normal (l)-vL[l]  [1] ♦Normal(O) ) ♦Normal (1))  ==0.0, 
vT  [0]  ♦Normal  (0)  +vT  [2]  ♦Normal  ( 1 )  +  (T  [2]  ♦Normal  (0  )  -T  [0]  ♦Normal  ( 1 )  )  ♦ 
((vL[0]  [0]  ♦Normal  (l)-vL[0]  [1]  ♦Normal  (0))  ♦Normal  (0)  + 

(vL[l]  [0]  ♦Normal  (l)-vL[l]  [1]  ♦Normal  (0))  ♦Normal  (l))==vPo[0]  && 
vT  [2]  ♦Normal  (0)  +vT  [1]  ♦Normal  ( 1)  +  (T  [1]  ♦Normal  (0)  -T  [2]  ♦Normal  ( 1) )  ♦ 
((vL[0]  [0]  ♦Normal (l)-vL[0]  [1]  ♦Normal  (0))^Nomal(0)  + 

(vL[l]  [0]^Normal(l)-vL[l]  [1] ♦Normal (0)) ♦Normal (l))==vPo[l]))  ; 
DiscExpr  elastSensit (vEquilEqns,grid,vU,vBC) ; 
vU=elastSensit . solve () ; 

} 
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Optimization  Appiication 

Square  Plate  with  Hole 


Optimize:  (p(p) 

9(P)=  ^  J  (aee(5)  -  O0e)^dy 


Figure  5:  Schematic  representation  of  optimization  problem. 


4.2.3  Application  of  PDESolve  to  Optimization  Problems  in  Elasticity 

As  an  application  of  optimization-based  design  in  linear  elasticity,  in  this  section  we  compute  the  optimal  shape 
of  the  cutout  in  a  thin  elastic  plate  subjected  to  constant  ratios  of  far  field  stresses  (see  Fig.  5).  Specifically, 
we  want  to  determine  the  optimal  shape  of  the  cutout  such  that,  for  the  given  far  field  stresses,  we  obtain  a 
uniform  distribution  of  tangential  stresses  along  the  boundary  of  the  cutout. 

The  particular  objective  function  to  be  minimized  in  this  problem  is  given  by 

(j>  =  ^  [  {(T00  -  a00fdl 

JdBc 

where  aee  is  the  tangential  stress  along  the  cutout,  and 

a00  =  —  f  a$$  dl  ,  Lc  =  f  dl 

I>c  JdBc  JdBc 

The  variation  of  </>  with  respect  to  the  design  variable  Pi  {5(f>  =  §^6Pi)  is  expressed  as 


/  ( 

JdBc 


((J00  —  (700){Sa00  5a00)dl 


+  ^  f  {(700  -  a00)Hdl  - -^SLc 
I'c  JdBb 


where 


J^c  JdBc 


c^ee 


i6a00dl  +  a00Sdl)-^5L^ 


6Lc  =  I  Sdl  ,  Sdl  =  {SLiiTi^  +  ~  iSLi2  +  S L2i)n,in2))dl 

JdBc 

The  problem  is  modeled  as  a  two-dimensional  plane  stress  problem.  By  symmetry,  only  one  quarter  of  the 
plate  is  analyzed.  The  material  properties  of  the  plate  are:  E  =  5.88  x  10^  MPa  (Young’s  modulus),  i/  =  0.3 
(Poisson  ratio).  The  geometry  and  boundary  conditions  of  the  problem  are  shown  in  Fig.  5.  The  applied  load 
is  taken  as  an  =  1.0  Mpa,  a22  =  0.75  Mpa.  The  analytical  solution  of  this  problem  [29]  shows  that  the 
optimum  shape  of  the  hole  is  elliptic,  with  the  semi-axes  a  and  b  given  the  ratio 


b 

a 


^22 

^11 


0.75 
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Figure  6:  Convergence  of  hole  geometry  to  optimum  shape.  One  design  parameter.  Linear  Elasticity. 


Figure  7:  VonMises  stress  for  initial  hole  geometry.  One  Figure  8:  VonMises  stress  for  optimized  hole  geometry, 
design  parameter.  Linear  Elasticity  One  design  parameter.  Linear  Elasticity 
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Figure  9:  Convergence  of  hole  geometry  to  optimum  shape.  Six  design  parameters.  Linear  Elasticity. 


Figure  10:  VonMises  stress  for  initial  hole  geometry.  Figure  11:  VonMises  stress  for  optimized  hole  geome- 
Six  design  parameters.  Linear  Elasticity  try.  Six  design  parameters.  Linear  Elasticity 
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The  PDESolve  code  for  solving  this  example  is  similar  to  the  one  presented  in  the  previous  section  and 
will  not  be  shown  here.  Two  cases  are  solved  with  this  code.  In  the  first  case,  the  shape  of  the  cutout  is 
assumed  to  be  elliptic  and  modeled  as  Xi  =  a  cos0,  X2  =  b  sin0.  The  geometry  of  the  problem  is  defined  by 
L  =  5m,  a  =  lm  and  b  is  the  design  parameter.  Numerical  results  are  presented  in  Figs.  6-8.  Fig.  6  shows 
the  different  geometries  of  the  hole  computed  during  the  iterative  solution  procedure.  The  starting  geometry 
was  a  circle  with  b  =  I  m  (</>  =  1.5  x  10“^),  and  the  final  (optimum)  geometry  is  an  ellipse  with  b  =  0.755 
m  (<^  =  4.7  X  10“®).  The  optimum  design  was  obtained  in  3  iterations.  Figures  7  and  8  shows  the  Von 
Mises  stress  distribution  throughout  the  plate.  Note  that  at  the  hole  boundary,  the  Von  Mises  and  tangential 
stresses  are  the  same.  Fig.  7  is  the  stress  distribution  prior  to  optimization  (6=1)  and  Fig.  8  is  the  stress 
distribution  after  optimization  (6  =  0.755).  It  is  observed  that  for  the  optimum  value  of  6  the  Von  Mises  stress 
is  distributed  uniformly  along  the  hole. 

For  the  second  case,  we  take  L  =  1  m,  a  =  1  m,  the  optimum  analytical  value  6  =  0.75  m,  and  start 
with  an  arbitrary  hole  geometry  represented  by  two  pieces  of  third  order  Bezier  curves  (14  control  points). 
The  selected  geometry  of  the  hole,  the  symmetry  conditions  of  the  plate  and  the  continuity  requirement 
imposed  at  the  joint  point  of  the  two  curves  leave  us  with  6  independent  control  points  which  are  taken 
as  the  shape  design  parameters.  The  initial  guess  for  these  parameters  defines  the  hole  geometry  presented 
in  Fig.  10.  This  figure  also  shows  the  distribution  of  Von  Mises  stress  throughout  the  plate  for  this  initial 
geometry.  Fig.  9  shows  the  different  geometries  of  the  hole  computed  during  the  iterative  solution  procedure. 
It  is  noted  that  the  optimum  shape  of  the  hole  (ellipse  or  close  to  an  ellipse)  is  obtained  after  8  iterations  (the 
objective  function  reduces  from  (j)  =  5.4  x  10“^  for  the  initial  geometry  to  ^  =  8.1  x  10“®  for  the  optimum 
geometry).  Fig.  11  shows  the  Von  Mises  stress  distribution  for  the  optimum  hole  shape.  It  is  observed  again 
that  for  this  optimum  shape,  the  stress  is  distributed  uniformly  along  the  hole. 

4.3  Design  Optimization  in  Small  Deformation  Plasticity 

4.3.1  Problem  Formulation 

For  the  general  formulation,  we  use  a  unified  elasto-viscoplastic  model  with  a  single  scalar  internal  variable 
(hardness)  to  characterize  the  plastic  deformation  of  a  material.  The  specific  (rate)  equations  are  given  by 


the  small  strain  version  of  the  constitutive  equations  (14)“(17), 

O'  =  C:(€-~€P)  (69) 

(70) 

f  =  fia,s)  (71) 

s  =  g{a,s)  (72) 


where  N  =  -0!^,  ^  =  y  |  II  o'*  II  tr'  =  <r  —  |ir(<r)I.  Here,  e  and  denote  the  total  strain  rate  and 
its  plastic  component,  respectively.  These  equations  together  with  eqs.(61)-(63)  (equilibrium  equations  and 
boundary  conditions)  define  the  direct  problem  for  small  deformation  plasticity.  The  numerical  solution  of 
these  nonlinear  equations  is  carried  out  using  an  incremental,  displacement-based  approach.  For  this  purpose, 
the  specified  time  interval  in  divided  into  a  certain  number  of  steps  A^,  with  a  typical  time  step  defined  by 
(tiv,  tyv+i)  where  i^v+i  =  +  A^.  Assuming  that  the  problem  has  been  solved  up  to  time  t^v,  the  incremental 

solution  procedure  is  focussed  on  computing  the  corresponding  field  variables  at  In+i-  As  before,  this  solution 
procedure  involves  iterations  at  both  the  equilibrium  and  the  constitutive  level.  For  the  equilibrium  iterations, 
we  write  eqs.(61)-(63)  at  In-vi  ^ 


V*(Tjv4-1  -f  f/v+1 

=  0  in 

V 

(73) 

=  tw+i 

on  Aa 

(74) 

UAT+i 

=  Ujv+l 

on  Au 

(75) 

For  algorithm  development,  it  is  convenient  to  express  the  stress  and  displacement  at  as:  (Tn^-i  = 
^{«)  -I- A^  and  u/v+i  =  -h  Au,  where  A<t  and  Au  are  the  incremental  stress  and  displacement,  respectively. 
Here,  the  i-superscripted  quantities  denotes  values  of  <t  and  u  at  the  i-th  iteration,  with  =  (Tn  and 
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=  uat.  Using  these  expressions  for  and  we  can  write  eqs(73)-(75)  as 

V-(A<t)  =  -f^+i  -  V  •  (76) 

Atr-n  =  In+i  — 

Au  =  un+i  — 

In  these  equation,  the  stress  increment  Aa  is  related  to  the  strain  increment  Ae  using  the  elasto- viscoplastic 
constitutive  matrix  £ 

Act  =  CiAe  (79) 

where 

C  =  2/i2:+(«-^/i)I®I-2M(»?-c)N(*^ON(‘)  (80) 

o 

Ae  =  i[V(Au)  +  (V(Au))’’]  (81) 

The  equilibrium  iterations  are  then  accomplished  using  eqs.(76)-(81)  to  obtain  the  incremental  displacement 
Au.  On  the  other  hand,  the  constitutive  iterations  are  carried  out  using  the  integrated  form  of  the  constitu¬ 
tive  equations,  eqs.(69)-(73),  to  update  the  stress  and  material  state  (hardness)  Siv+i-  Details  of  the 

constitutive  integration  procedure  are  presented  in  [22] 

For  the  sensitivity  problem,  the  governing  equations  are  given  by 


in  V 

(82) 

(^o'Tv+i  “■  (JL:n  (8)  n)tjv+i  = 

on  Aa^ 

(83) 

on  Au 

(84) 

together  with  the  sensitivity  equations  for  the  constitutive  model.  These  last  equations  relates  the  stress 
variation  Scr^^i  to  the  strain  variation  5Ae  and  are  derived  from  the  integrated  form  of  the  eqs.(69)-(73). 


Omitting  details  of  this  derivation,  these  equations  are 

(5<7iv+i  =  £jv+i'<5Ae  +  *AiVH-i  (85) 

(5siv+i  =  —  (6i(5a*  +  aiJsiv)  (86) 

h 

where  Sctn  and  Ssn  represent  the  variation  of  the  Cauchy  stress  and  the  hardness  at  time  tjv?  respectively, 
and 

An+1  =  vil-v  g  I  c)N7v+i  ®  Njv+1  (87) 

6aN  = 

SAe  =  i[V(5Au+(V<5Au)^-(V(Au)-(5L  +  (V(Au)-<5L)^)]  (89) 


Equations  (82) -(89)  are  solved  for  (5Au  at  each  time  step  once  the  solution  for  the  direct  problem  had  been 
obtained.  Note  that  the  sensitivity  equations  are  linear  in  SAu 


4,3.2  Implementation  of  2-D  Small  Deformation  Plasticity  in  PDESolve 

The  previous  equations  have  been  simplified  to  2-Dimensional  plane  stress  and  plane  strain  cases,  and  imple¬ 
mented  using  PDESolve.  The  module  for  the  solution  of  the  direct  plane  strain  problem  is  presented  below. 
The  prescribed  boundary  conditions  shown  in  this  program  correspond  to  those  of  a  thick- wall  cylinder  under 
a  prescribed  velocity  at  its  inner  surface.  The  analytic  solution  for  the  radial  stress  (Xrr  is  given  in  [26]  for 
a’  rigid,  rate-independent  plastic  material.  This  solution  and  its  corresponding  sensitivity  Sarr  to  the  inner 
radius  a  are  compared  here  with  the  values  predicted  by  the  PDESolve  program.  The  particular  viscoplastic 
constitutive  model  used  in  this  analysis  is  a  power  law  of  the  form 

a  =  <To(l  +  -)"(1  + 

€o  Co 


(90) 
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where  a  is  the  equivalent  (Von  Mises)  stress,  is  the  equivalent  plastic  strain,  P  is  the  equivalent  plastic 
strain  rate,  and  cro,CO)€o,n,m  are  material  parameters.  In  this  example,  we  use:  ao  =  50  MPa,  eo  =  0.002, 
60  =  0.0025*“^,  n  =  0,  m  =  0.  The  elastic  constants  are  E  =  25  x  10^  MPa,  u  =  0.3.  The  problem  is  solved 
incrementally,  using  24  increments  with  a  constant  time  step  of  At  =  0.5s.  The  prescribed  velocity  at  the 
inner  surface  is  ^  =  1.0.  Figure  12  presents  the  nondimensional  radial  displacement  ^  versus  pressure 
ratio  ^  (p  =  ““^rr)  for  a  material  point  located  at  the  inner  surface  of  the  cylinder.  The  numerical  solution 
was  computed  using  both  second  (D2)  and  fourth  (D4)  order  of  accuracy  in  the  representation  of  the  finite 
difference  operators  =  1,2.  The  results  are  in  good  agreement  with  the  analytical  solution.  Figure  13 

shows  ^  versus  for  the  same  material  point.  As  before,  the  symbols  S2  and  S4  indicate  the  order  of 
accuracy  of  the  finite  difference  stencil  used  to  solve  the  sensitivity  problem.  It  is  important  to  note  the  effect 
of  the  finite  difference  representation  on  the  computed  sensitivities,  in  particular  in  the  transition  region  from 
elastic  to  viscoplastic  behavior.  As  deformation  proceeds,  the  calculated  sensitivities  have  the  same  trend  as 
the  analytical  ones,  obtaining  for  this  problem  a  better  match  with  the  analytical  solution  using  the  fourth 
order  precision  of  the  finite  difference  representation. 

void  Small_Plasticity.Direct () 

f 

//  Domain  and  Finite  Difference  Mesh 

SpatialDomain  computationalDomain  =  SpatialDomain(0, l)*SpatialDomain(0, 1) ; 

Function  fwdMap (2, Tensor (2,1) ,fwdMapFunc .Function (2, Tensor (2, 2) ,DfwdMapFunc)  ) ; 

Function  invMap(2,Tensor(2,l) ,invMapFunc) ; 

SpatialDomain  physicalDomain (computationalDomain,  Map(fwdMap,  invMap)); 

FDMesh  grid  =  FDMesh (physicalDomain,  FDMeshSpec (nX)  ♦  FDMeshSpec (nY) ) ; 


//  Differential  Operators 
DiffOp  Dx (grid, 1,0, 4); 
DiffOp  Dy(grid,l,l,4) ; 


//  Function  Definitions 

Function  U(2,Tensor(2,l)) ; 

Function  dU(2,Tensor(2,l)) ; 

Function  T(2,Tensor(4,l)) ; 

Function  EE (2, Tensor (4, 1) ) ; 

Function  EP(2,Tensor(4,l)) ; 

Function  EQP(2,Tensor (1, 1)) ; 

Function  Tt(2,Tensor(4,l)) ; 

Function  EEt(2,Tensor(4,l)) ; 

Function  EPt (2,Tensor (4, 1)) ; 

Function  EQPt (2, Tensor (1 , 1)) ; 

Function  J4x4(2,Tensor(4,2)) ; 

Function  F(2,Tensor(2, 1) ,  bodyForceFunc) ; 

Function  velBar (2, Tensor (2,1) ,BCVelocFuncMag) ; 

Function  UBar(2,Tensor(2, 1) ,zeroVectorFunc) ; 

Function  strikeVect4To3 (2, Tensor (3,1) ,strikeFuncVect4To3) ; 
Function  strikeMatr4To3 (2, Tensor (3,2) ,strikeFuncMatr4To3) ; 


/♦displacement*/ 

/♦accumulated  incremental  displacement*/ 
/♦Cauchy  stress*/ 

/♦elastic  strain*/ 

/♦plastic  strain*/ 

/♦equivalent  plastic  strain*/ 

/♦trial  Cauchy  stress*/ 

/♦trial  elastic  strain*/ 

/♦trial  plastic  strain*/ 

/♦trial  equivalent  plastic  strain*/ 
/♦Consistent  Tangent  Operator*/ 

/♦body  force*/ 

/♦prescribed  velocity*/ 

/♦accumulated  prescribed  displacement*/ 
/♦vector  4x1  ->  vectpr  3x1*/ 
/♦matrix  4x4  ->  matrix  3x3*/ 


//  Problem  Solution 
U  =  zeroVec; 

T  =  zeroTen; 

EE  =  zeroTen; 

EP  =  zeroTen; 

EQP  =  zeroSca; 

for  (incr  =  0;  incr  <  nincrs;  ++incr)  { 
dU  =  zeroVec; 

UBar  =  UBar  +  dTime  ♦  velBar; 
int  itersEquil  =  0; 
do  { 

if  (itersEquil  >  maxItersEquil)  throwRunTimeErr or ("Equilibrium  did  not  converge"); 
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//  Constitutive  Model  Integration  and  Consistent  Tangent  Operator 

Tt  =  T; 

EEt  ^  EE I 
EPt  =  EP; 

EQPt  =  EQP; 

Expr  diffOperDirect  =  List (List (Dx,  0.0,  Dy) ,  List (0,0,  Dy,  Dx)); 

Function  dE3  =  transpose (diffOperDirect) ♦dU; 

Function  dE4  =  List(dE3[0]  ,  dE3[l],  0.0,  dE3[2]); 

Function  listlnMat  =  List(dE4,  Tt,  EEt,  EPt,  EQPt); 

Function  listOutMat  =  MaterialConstIntegOp (listlnMat) ; 

Tt  =  listOutMat [0] ; 

EEt  =  listOutMat [1] ; 

EPt  =  listOutMat [2] ; 

EQPt  =  listOutMat [3] ; 

J4x4  =  listOutMat [4] ; 

Function  J3x3  =  strikeMatr4To3 ( J4x4) ; 

Function  Ttt  =  strikeVect4To3(Tt) ; 

//  Compute  the  displacement  increment  ddU 

Fxmction  ddU(2,  Tensor (2,  D); 
ddU . regardAsUnknown  0 ; 

Expr  ddE  =  transpose (diffOperDirect )*ddU; 

Expr  dT  =  J3x3  ♦  ddE; 

Expr  equilEqns  =  dif f OperDirect* (dT  +  Ttt)  +  F; 

BC  be  =  BC (List ((dT[0]+Ttt[0]) ♦Normal (0)  +  (dT[2]+Ttt [2] )+Normal(l)  ==  0.0, 
U[l]+ddU[l]  ==  0.0), 

List(U[0]  +  ddU[0]  ==  0.0, 

(dT[2]+Ttt[2])^Normal(0)  +  (dT[l]+Ttt  [1]  )  ♦Normal  (1)  ==  0.0)) 

♦  BC(List(U[0]  +  ddU[0]  ==  UBar[0], 

U[l]  +  ddU[l]  ==  UBar[l]), 

List  ((dT[0]+Ttt[0])  ♦Normal  (0)  +  (dT[2]+Ttt  [2]  )^Normal(l)  ==  0.0, 
(dT[2]+Ttt[23)^Normal(0)  +  (dT[l]+Ttt  [1]  )  ♦Normal  (1)  ===  0.0)); 
DisePDO  plastDirect (equilEqns ,  grid,  ddU,  be); 
ddU  =  plastDirect . solve () ; 

//  Update  U  and  dU;  Check  Convergence  in  ddU 

dU  =  dU  +  ddU; 

U  =  U  +  ddU; 
itersEquil++ ; 
errl  =  12Norm(ddU) ; 
err2  =  12Norm(U) ; 

> 

while  ( (errl  /  err2)  >  tolerEquil) ; 

//  Update  Veiriables 
T  =  Tt; 

EE  =  EEt; 

EP  =  EPt; 

EQP  =  EQPt; 


4,3.3  Application  of  PDESolve  to  Optimization  Problems  in  Small  Deformation  Plasticity 

In  this  section,  we  solve  the  shape  optimization  of  elliptic  cutouts  in  thin  elasto-viscoplastic  plates  subjected 
to  constant  ratios  of  far  field  stresses  (see  Fig.  5).  Again,  we  want  to  compute  the  optimum  shape  of  the  cutout 
for  a  uniform  tangential  stress  distribution  along  it.  Hence,  the  objective  function  and  its  sensitivity  are  the 
same  as  the  ones  presented  before.  In  this  example,  the  deformation  of  the  plate,  modeled  as  a  plane  stress 
problem,  is  driven  by  the  time  dependent  loading:  (Txxit)  =  16.0  +  8^  MPa,  (Tyy(t)  =  0.75axx(i)  MPa,  where  t 
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Figure  12:  Radial  stress  (arr  =  -p)  at  r  =  a  of  thick  Figure  13:  Sensitivity  Sarr  at  r  =  a  of  thick  wall  cylin- 
wall  cylinder.  Small  Deformation  Plasticity  der.  Small  Deformation  Plasticity 


Figure  14:  Convergence  of  hole  geometry  to  optimum  shape.  Small  Deformation  Plasticity. 


Figure  15:  VonMises  stress  for  initial  hole  geometry.  Figure  16:  VonMises  stress  for  optimized  hole  geome- 
Small  Deformation  Plasticity  try.  Small  Deformation  Plasticity 
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is  time.  We  assume  an  elliptic  geometry  of  the  hole  (rri  =  acos0,  X2  =  bsind)  and  take  L  =  5  m,  a  =  1  m,  with 
b  being  the  design  parameter.  The  material  parameters  used  are  (Tq  =  50  MPa,  €o  ==  0.002,  eo  =  0.0025"^^, 
n  =  0.1,  m  =  0.1,  £?  =  25  X  10^,  1/  =  0.3.  The  plasticity  solution  of  the  direct  and  sensitivity  problems  is 
carried  out  incrementally  for  a  prescribed  time  interval.  Specifically,  the  problem  is  solved  using  25  increments 
with  a  constant  time  step  of  At  =  0.16s  to  give  a  total  time  of  4s.  We  use  a  second  order  accuracy  of  the 
finite  difference  representation  for  solving  both  the  direct  and  the  sensitivity  problems.  Numerical  results 
are  presented  in  Figs.  14  -  16.  Fig.  14  shows  the  different  hole  geometries  computed  during  the  optimization 
process.  The  starting  geometry  was  a  circle  with  6  =  1  m  (<^  =  59.7),  and  the  final  (optimum)  geometry  is 
an  ellipse  with  b  =  0.667  m  (^  =  7.6  x  10”^).  The  optimum  elasto-viscoplastic  design  is  obtained  after  4 
iterations.  Note  here  that  the  stress  distribution  along  the  cutout  for  the  elasticity  analytical  solution  b  =  0.75 
m  is  not  longer  the  optimal  when  inelastic  deformation  occurs.  Hence,  the  optimum  shape  of  the  cutout  for 
the  new  stress  pattern  is  given  by  6  =  0.667  m.  Figs.  15  and  16  shows  the  Von  Mises  stress  distribution 
throughout  the  plate.  Fig.  15  is  the  stress  distribution  prior  to  optimization  (6  =  1)  and  Fig.  16  is  the  stress 
distribution  after  optimization  (6  =  0.667).  We  observe  that  for  the  optimum  value  of  b  the  Von  Mises  stress 
distribution  is  uniform  along  the  hole. 

5  PDESolve  AApplied  to  Shape  Optimization  of  One-Step  Open 
Die  Forging 

5.1  Implementation  of  2-D  Large  Deformation  Plasticity  in  PDESolve 

Before  applying  the  general  optimization-based  design  methodology  to  the  forging  problem,  in  this  section  we 
validate  the  2-D  plane  strain  implementation  of  large  deformation  plasticity  using  PDESolve.  The  particular 
module  for  solving  the  large-deformation  direct  problem  is  presented  below.  With  this  code,  we  solve  the  direct 
and  sensitivities  problems  of  the  large  deformation  of  a  thick-wall  cylinder  subjected  to  an  applied  velocity  at 
its  inner  surface.  As  before,  the  analytical  solution,  given  in  [26],  is  used  here  to  compare  with  the  prediction 
of  PDESolve.  We  use  the  viscoplastic  model  given  by  eq.(90)  with  the  material  parameters:  cro  =  50  MPa, 
€o  =  0.002,  eo  =  0.002s“^  n  =  0,  m  =  0.  Also,  E  =  25  x  10^  MPa  and  u  =  0.3.  The  problem  is  solved 
incrementally,  using  25  increments  with  a  constant  time  step  of  At  =  20.0s  to  give  a  total  time  of  500s.  The 
prescribed  velocity  at  the  inner  surface  is  again  ^  =  1.0.  The  plots  ^  vs  ^  and  ^  vs  for  a  material 
point  at  the  inner  surface  of  the  cylinder  are  shown  in  Figs.  17  and  18,  respectively.  We  observe  that 
the  numerical  solution  of  the  direct  problem  is  in  good  agreement  with  the  analytical  one.  The  computed 
sensitivity  coeficients,  on  the  other  hand,  have  the  same  trends  as  the  analytical  ones,  with  a  maximum  error 
of  approximately  25%. 

void  Large_Plasticity_Direct() 

f 

//  Domain  and  Finite  Difference  Mesh 

SpatialDomain  computationalDomain  =  SpatialDomain(0,l)*SpatialDomain(0, 1) ; 

Function  fwdMap (2, Tensor (2,1) ,fwdMapFunc, Function (2, Tensor (2,2) ,DfwdMapFunc) ) ; 

Function  invMap(2,Tensor (2,1) ,invMapFunc) ; 

SpatialDomain  undeformedDomain( computationalDomain,  Map (fwdMap, invMap)) ; 

FDMesh  undeformedMesh  =  FDMesh(undeformedDomain,  FDMeshSpec(nX)*FDMeshSpec(nY)) ; 

FDDiscretize  undeformedBody (undeformedMesh) ;  /*initial  configuration,  B_0*/ 

Function  X (2, Tensor (2,1) ,coordFunc) ; 

Function  Xini  =  undeformedBody .discretize (X) ;  /*initial  coordinates,  x„0*/ 

Function  Xref  =  undeformedBody. discretize (X) ;  /*reference  coordinates,  x_N*/ 

Function  Xcur(2,Tensor(2,l)) ;  /^current  coordinates,  x_(N+l)*/ 

//  Differential  Operators 
DiffOp  Dx(l,  0,  4); 

'  DiffOp  Dy(l,  1,  4); 

Expr  Del  =  List(Dx,  Dy) ; 

//  Function  Definitions 

Function  U(2,Tensor(2,l)) ; 

Function  Ut(2,Tensor(2,l)) ; 


/♦displacement*/ 
/♦trial  displacement*/ 
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Function  dU(2,Tensor(2, 1) ) ; 

Function  T(2,Tensor(4,l)) ; 

Function  EE (2, Tensor (4, 1) ) ; 

Function  EP(2,Tensor(4, 1) ) ; 

Function  EQP(2,Tensor(l,l)) ; 

Function  DEQPt (2, Tensor (1 , 1)) ; 

Function  Tt(2,Tensor(4,l) ) ; 

Function  EEt(2,Tensor(4,l)) ; 

Function  EPt(2,Tensor(4,l)) ; 

Function  EQPt (2, Tensor (1,1)) ; 

Function  J4x4(2,Tensor(4,2)) ; 

Function  FN(2,Tensor(2,2)) ; 

Function  FNl(2,Tensor(2,2)) ; 

Function  DF(2,Tensor (2,2)) ; 

Function  dF(2,Tensor(4,l)) ; 

Function  J5x5(2,Tensor(5,2)) ; 

Function  B (2, Tensor (2, 1) ,bodyForceFunc) ; 

Function  velBeo: (2, Tensor (2,1) ,BCVelocFuncMag) ; 

Function  UBar (2, Tensor (2,1) ,zeroVectorFunc) ; 

//  Problem  Solution 

U  =  undeformedBody.discretize(zeroVec) ; 

T  =  undeformedBody.discretize(zeroTen) ; 

EE  =  undeformedBody. discretize (zeroTen) ; 

EP  =  undeformedBody, discretize (zeroTen) ; 

EQP  =  undef ormedBody, discretize (zeroSca) ; 
dTimeN  =  dt[0];  Z+This  module  includes  an  automatic  time  stepping  scheme*/ 
time  =  0.0;  /♦Not  used  in  this  problem*/ 

incr  "  0; 
counter  =  0; 
do  { 

incr  =  incr  +  1; 

//  Define  reference  configuration  and  setup  equilibrium  iterations 
SpatialDomain  ref erenceDomain (computat ionalDomain ,  Map (Xref ) ) ; 

FDMesh  referenceMesh  =  FDMesh (ref erenceDomain,  FDMeshSpec (nX)  *  FDMeshSpec(nY)) ; 
FDDiscretize  ref erenceBody (referenceMesh) ; 

FN  =  transpose (outerProduct (Del,  Xref)); 

Tn  =  undeformedBody .discretize (T) ; 
redoEquilSol :  ; 
convergeEquilSol  =  true; 
dTime  =  dTimeN; 

UBarTr  =  UBar  +  dTime  *  velBar; 

dU  =  referenceBody. discretize (zeroVec) ; 

Xcur  =  undeformedBody .discr etize( Xref ) ; 

Ut  =  undeformedBody. discretize (U) ; 
timeTr  =  time  +  dTime; 
int  itersEquil  =  0; 
do  { 

if  (itersEquil  >  maxItersEquil) 

{  convergeEquilSol  =  false; 

cout  «  "  convergeEquilSol  =  FALSE  for  dTime  =  ’*  «  dTime  «  endl; 
cout  «  "  ...  will  try  dTimeN  =  0 . 5*dTime  "  «  endl ; 
goto  checkTimeStep;  } 

//  Constitutive  Model  Integration  and  Consistent  Tangent  Operator 

Tt  =  referenceBody .discretize (T) ; 

EEt  =  referenceBody .discretize (EE) ; 

EPt  =  referenceBody .discretize (EP) ; 

EQPt  =  referenceBody .discretize (EQP) ; 

FNl  =  transpose (outerProduct (Del,  Xcur)); 


/♦accumulated  incremental  displacement*/ 
/♦Cauchy  stress*/ 

/♦elastic  strain*/ 

/♦plastic  strain*/ 

/♦equivalent  plastic  strain*/ 
/♦incremental  EQP*/ 

/♦trial  Cauchy  stress*/ 

/♦trial  elastic  strain*/ 

/♦trial  plastic  strain*/ 

/♦trial  equivalent  plastic  strain*/ 
/♦consisten  tangent  operator*/ 
/♦deformation  gradient  F_N*/ 
/♦deformation  gradient  F_(N+l)=F„u*/ 
/♦relative  deformation  gradient  F*/ 

/♦  F  in  vector  form*/ 

/♦consistent  tangent  operator*/ 

/♦body  force,  B=0  in  this  example*/ 
/♦prescribed  velocity*/ 

/♦accumulated  prescribed  displacement*/ 
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DF  =  FNl  ♦  inverseMatxOp(FN) ; 

dF  =  List(DF[0]  [0]  ,DF[0]  [1]  .DF[1]  [0]  ,DF[1]  [1])  ; 

Function  listlnMat  =  List (dF,Tt ,EEt,EPt ,EQPt) ; 

Function  listOutMat  =  MaterialConstIntegOp (listlnMat) ; 

It  =  listOutMat [0] ; 

EEt  =  listOutMat [1] ; 

EPt  =  listOutMat [2] ; 

EQPt  =  listOutMat [3] ; 

J5x5  =  listOutMat [4] ; 
listPlast  =  listOutMat [6] ; 

Function  J4x4  =  List  (List  (J5x5[0]  [0]  ,  J5x5[0]  [1]  ,  J5x5[0]  [2]  ,J5x5[0]  [3]  ) , 

List(J5x5[l]  [0]  ,  J5x5[l]  [1]  .J5x5[l]  [2]  ,J5x5Cl]  [3])  . 

List(J5x5[2]  [0]  .  J5x5[2]  [1]  ,J5x5[2]  [2]  ,J5x5[2]  [3])  . 

List  (J5x5  [3]  [0]  ,  J5x5  [3]  [1]  ,  J5x5  [3]  [2]  ,  J5x5[3]  [3]  )  )  ; 

Function  matxTl  =  List(List(Tt[0] ,0.0,0.0,Tt[0]) , 

List(Tt[l]  ,0.0,0.0,Tt[l]) , 

List(Tt[2],0.0,0.0,Tt[2]), 

List(Tt[3],0.0,0.0,Tt[3])); 

Function  matxT2  =  List(  List (Tt [0] ,0.0, Tt[l] ,0.0) , 

List  (0 . 0 ,  (0 . 5*  (Tt  [0]  -Tt  [3]  )  )  ,  (0 . 5*  (Tt  [0]  +Tt  [3] ) )  ,  Tt  [1]  )  , 
List  (Tt  [2]  .  (0 . 5*  (Tt  [0]  +Tt  [3]  )  )  ,  (0 . 5*  (-Tt  [0]  +Tt  [3]  ) ) ,  0 . 0)  , 
List(0.0,Tt[2] .0.0,Tt [3])); 

Function  matxX  =  List(  List(l .0,0. 0,0. 0,0.0) , 

List (0.0, 0.5, 0.5, 0.0) , 

List (0.0, 0.5, 0.5, 0,0) , 

List(0. 0,0. 0,0. 0,1.0)  ); 

Function  Ttt  =  List  (Tt[0]  ,Tt[l]  ,Tt  [2]  ,Tt[3])  ; 

//  Define  current  condifuration  and  compute  the  displacement  increment  ddU 

SpatialDomain  cur rent Domain ( comput at ionalDomain,  Map(Xcur)); 

FDMesh  currentMesh  =  FDMesh( current Domain,  FDMeshSpec (nX)  *  FDMeshSpec (nY) ) ; 
FDDiscretize  cvirrentBody( currentMesh)  ; 

Function  ddU(2,Tensor (2,1)) ; 
ddU . regardAsUnknown  0 ; 

Expr  diffOperDirect  =  List (List (Dx,Dy, 0.0, 0.0) , List (0.0, 0.0, Dx,Dy) ) ; 

Expr  ddF  =  transpose (diffOperDirect)  ♦  ddU; 

Expr  dP  =  (J4x4*matxX  -  matxT2  +  matxTl)  ♦  ddF;  /*!  Piola-Kircchhof  stress*/ 

Expr  equilEqns  =  diffOperDirect* (dP  +  Ttt); 

BC  be  =  BC(List(Ut[0]  +  ddU[0]  ==  UBarTr[0], 

Ut[l]  +  ddUCl]  ==  UBarTrEl]), 

List((dP[0]  +  Ttt[0])*Normal(0)  +  (dP[l]  +  Ttt [1] )*Normal(l)  ==  0.0, 

(dP[2]  +  Ttt[2])*Normal(0)  +  (dP[3]  +  Ttt  [3]  )*Normal(l)  ==  0.0)) 

*  BC(List((dP[0]  +  Ttt[0])*Normal(0)  +  (dP[l]  +  Ttt [1] )*Normal(l)  ==  0.0, 

Ut[l]  +  ddU[l]  ==  0.0), 

List(Ut[0]  +  ddU[0]  ==  0.0, 

(dP [2] +Ttt [2] ) *Normal (0)  +  (dP [3] +Ttt [3] ) *Normal (1)  ==  0.0)); 
DiscExpr  plastDirect (currentMesh,  equilEqns,  be,  ddU) ; 
ddU  =  plastDirect .solve 0 ; 

//  Update  Ut  and  dU;  Check  Convergence  in  ddU 

Xcur  =  undef ormedBody. discretize (Xcur  +  ddU) ; 
dU  =  referenceBody. discretize (dU  +  ddU) ; 

Ut  =  undef ormedBody .discretize (Ut  +  ddU) ; 

itersEquil++; 

errl  =  12Norm(ddU) ; 

err2  =  12Norm(Ut) ; 

} 

while  ( (errl  /  err2)  >  tolerEquil) ; 


//  Compute  new  or  next  time  step  (For  Automatic  Time  Stepping  Scheme) 
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checkTimeStep:  ; 

DEQPt  =  referenceBody. discretize (listPlast) ; 
if  (autoTime  ==  0) 

{  dTimeN  =  dt [incr] ; 
rejectEquilSol  =  false; 
if  (IconvergeEquilSol) 

throwRunTimeErrorC'No  autoTime  was  prescribed  ;  exit  program”);} 
else  if  (autoTime  ==  1) 

{  timeControlFuncCconvergeEquilSol,  rejectEquilSol,  dTime, 
dTimeN,  counter,  DEQPt) ;  } 

if  (rejectEquilSol) 

{  cout  «  ”  WARNING;  Redoing  Equilibrium  Solution"  «  endl; 
goto  redoEquilSol ;  } 

//  Update  Body’s  Configuration  and  Variables 
Xref  =  undeformedBody. discretize (Xcur) ; 

SpatialDomain  cur rent Domain ( comput at ionalDomain,  Map (Xref)); 

FDMesh  currentMesh  =  FDMesh(currentDomain,  FDMeshSpec(nX)  *  FDMeshSpec (nY) ) ; 
FDDiscretize  currentBody( currentMesh) ; 
time  =  timeTr; 

UBar  =  UBaxTr; 

U  =  undeformedBody .discretize (Ut) ; 

T  =  currentBody.discretize(Tt) ; 

EE  =  currentBody.discretize(EEt) ; 

EP  =  currentBody.discretize(EPt) ; 

EQP  =  currentBody. discretize (EQPt) ; 

} 


5.2  Shape  Optimization  of  One-Step  Open  Die  Forging 

In  this  section,  we  use  the  general  methodology  developed  for  optimization-based  design  of  multi-step  forging 
to  solve  the  preform  shape  design  of  one-step  forging  (upsetting)  where  a  rigid  die  of  fixed  geometry  is  used. 
For  this  particular  case,  the  methodology  is  simplified  to  the  following  steps: 

1.  Guess  the  initial  value  of  the  design  parameters  ft. 

2.  Do  i  =  1,  #  steps  (incremental  analysis): 

(a)  Solve  direct  problem  for  the  field  variables  (u,  P).  Use  the  integrated  form  of  the  constitutive 
equations,  eqs.(20)-(23),  and  the  incremental  equilibrium  equations,  eqs.(27)-(29). 

(b)  Solve  n  sensitivity  problems  for  the  variations  (Ju,  SP),  where  n  is  the  number  of  design  parameters. 
Use  eqs.(35)‘-(37)  together  with  eqs.(51)-(52). 

3.  Solve  optimization  problem  to  compute  new  values  of  /3.  Here,  the  objective  function  ^  (see  below)  is 
minimized  using  gradient  information  (S^). 

4.  Check  convergence  of  the  objective  function: 

If  II  $  II <  Tolerance  =>  exit 

If  11  #  11  >  Tolerance  =>  go  to  step  2  with  the  new  values  of  /3. 

The  specific  problem  to  be  solved  is  the  design  of  the  preform  shape  such  that  after  the  forging  (upsetting) 
process,  the  final  product  has  a  minimum  barreling  (see  fig.  19).  Modeling  the  problem  as  a  two-dimensional 
plane  strain  problem  (only  one  quarter  of  the  preform  need  to  be  considered),  we  can  write  the  objective 
function  as  [2] 

$(/3)=  f  (xi-ro)"dr  (91) 

Jr, 

where  tq  is  a  prescribed  value,  F/  is  the  free  surface  of  the  deformed  workpiece  after  upsetting,  and  the 
coordinate  xi  =  x-ei  locates  a  material  point  on  this  free  surface.  Here,  x  =  xq  H-  u,  and  ei  is  a  unit  vector 
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Figure  17:  Radial  stress  (orr  =  —p)  at  r  =  a  of  thick  Figure  18:  Sensitivity  Sarr  at  r  =  a  of  thick  wall  cylin- 
wall  cylinder.  Large  Deformation  Plasticity  der.  Large  Deformation  Plasticity 


Figure  20:  Convergence  of  preform  free  surface  geometry  to  optimum  shape.  Four  design  parameters.  One-step 
forging. 


BEAM  Technologies,  Inc.  PROPRIETARY 


31 


along  the  1-direction,  The  variation  of  $  is  computed  as 

6^  —  f  [2{xi  —  ro)5xi  +  {xi  —  ro)^(<r((5L)  —  <5L:n  0  n)]dr 
hs 


(92) 


where 


6L  =  6Ft^F:^^  =SF^-FZ^ ,  SLo  =  6Fr^F^^  (93) 

Sx  =  ^xo  +  ^u  ,  Jxi  =  Jx-ei  ,  5xoi  =  <5xoi*ei  (94) 

In  these  expressions,  the  field  variables  (u,  F^)  are  computed  from  the  solution  of  the  direct  problem  while 
the  variations  (Ju,  5F„)  axe  obtained  from  the  sensitivity  analysis  of  the  problem.  As  shown  in  fig.  19,  the 
preform  geometry  is  mapped  to  the  2-D  unit  square  using  the  bilinear  transfinite  mapping 

Xox  =  mF{V2) 

Xo2  =  ViG{V2)  +  mil  -  m)L 

where  the  functions  F{7]2)  and  G{rj2)  define  the  geometry  of  the  free  boundary  of  the  preform  by  a  third  order 
Bezier  curve,  i.e., 

^(772)  =  (1  ^  V2)^^0  +  37/2(1  —  +  37/1(1  —  7/2)02  +  7/2 03 

Gim)  =  (1  -  +  37/2(1  -  772)^61  +  37/2(1  -  7/2)62  +  r]lh 

where  ai,bi,i  =  0,3  are  the  control  points  of  the  Bezier  curve  with  F(0)  =  6,  F{1)  =  c,  G(0)  =  0,  G(l)  =  L 
(see  Fig.  19).  This  mapping  is  used  to  compute  the  variation  of  the  design  velocity  gradient  <5Lo,  which  for  the 
case  of  preform  shape  design,  drives  the  sensitivity  problem.  For  the  numerical  simulations,  we  will  assume 
sticking  friction  at  the  workpiece-die  interface.  Also,  we  will  use  the  viscoplastic  relation  defined  by  eq.(90) 
with  the  material  parameters  ao  =  50  MPa,  eo  =  0.002,  eo  =  0.002s“^,  n  =  0.1,  m  =  0.1,  E  =  25x  10^  MPa, 
1/  =  0.3.  The  workpiece  is  deformed  up  to  20%  upsetting  with  a  compressive  velocity  of  12  x  10“^  mm/s  using 
25  increments  with  a  time  step  size  of  At  =  1.0s  (nominal  strain  rate  of  0.009s“^) 

The  geometry  of  the  preform  and  the  symmetry  conditions  give  ao  =  ai  =  6,  as  =  c,  60  =  0  and  63  =  L.  For 
this  example,  we  fixed  the  values  c  =  1  mm,  L  =  1.5  mm,  ro  =  1.0  mm  and  take  ao(=  0i),a2i=  ^2)j6i(=  (3s) 
and  62 (=  /?4)  SIS  the  four  design  parameters  ft.  The  initial  guess  for  these  parameters,  given  by:  ft  =  1.0, 
ft  =  1.0,  ft  =  0.5,  and  ft  =  1.0,  defines  the  initial  preform  geometry  shown  in  Fig.  21.  Figure  22  shows 
the  distribution  of  equivalent  plastic  strains  throughout  the  deformed  workpiece  after  20%  upsetting.  We 
clearly  see  the  barreling  of  the  workpiece  due  to  friction.  For  this  deformed  geometry,  the  objective  function 
has  the  value  </>  =  8.3  x  10“^.  Figure  20  presents  the  different  shapes  of  the  preform  free  surface  computed 
during  the  optimization  process.  The  ’’optimum”  shape  is  obtained  after  5  iterations.  The  optimum  design 
parameters  are  ft  =  ft  =  0.673,  ft  =  0.496,  ft  =  0.886,  ft  =  0.979  giving  a  value  for  the  objective  function  of 
^  =  4.3  X  lO”'*.  The  initial  preform  geometry  and  the  corresponding  deformed  workpiece  are  shown  in  Figs.  23 
and  24.  Although  the  deformed  workpiece  does  not  have  a  perfectly  vertical  free  surface,  it  is  observed  that 
the  geometry  of  the  forged  product  is  close  to  the  one  specified.  We  note  here  that  this  solution  can  be 
improved  by  using  a  better  representation  of  the  free  surface  as  well  as  by  mesh  refinement. 

6  A  Roadmap  for  Constitutive  Models  to  be  Included  in  the  Frame¬ 
work 

The  framework  being  developed  in  Phase  I  for  optimization-based  forging  design  uses  a  hypo-elastic  constitu¬ 
tive  model  [22]  for  the  material  representation  and  simple  contact /friction  conditions  (such  as  sticking  friction) 
at  the  workpiece-die  interface.  Such  simplified  framework  has  been  coded  using  our  computational  substrate 
PDESolve  to  solve  preform  design  problems  for  one-step  isothermal  forging.  For  Phase  II,  this  framework 
will  be  extended  to  incorporate  constitutive  laws  that  appropriately  represent  the  thermo-elasto-viscoplastic- 
damage  material  behavior  of  typical  metallic  alloys  of  interest  to  the  Air  Force.  A  number  of  material  models 
together  with  their  numerical  time  integration  procedure  will  be  investigated.  Appropriate  models  for  repre¬ 
senting  the  contact /friction  conditions  at  the  workpiece-die  interface  as  well  as  the  corresponding  numerical 


Figure  23:  Preform  geometry  for  optimum  values  of  Figure  24:  Deformed  workpiece  for  optimum  preform 
design  parameters.  One-step  forging.  geometry  after  20%  upsetting.  One-step  forging. 
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algorithms  will  also  be  addressed.  This  extended  constitutive  description  will  give  a  more  accurate  description 
of  the  physics  that  controls  the  final  product  in  a  forging  process. 

One  such  aspect  of  the  physics  is  the  formation  of  voids  or  cavities  in  regions  where  high  tensile  stresses  are 
produced.  This  phenomena  leads  to  degradation  of  mechanical  properties  of  the  forged  product,  and  hence, 
to  its  premature  failure.  We  will  incorporate  damage  models  into  the  framework  to  predict  zones  where  such 
conditions  may  occur  and  use  this  information  in  the  optimization  process  to  minimize  the  risk  of  void  forma¬ 
tion.  Specific  models  based  on  phenomenological  continuum  approaches  that  account  for  isotropic/anisotropic 
damage  descriptions  will  be  investigated  [13,  35]. 

Another  aspect  of  the  physics  is  the  nonuniform  distribution  of  plastic  strain  throughout  the  workpiece 
imparted  during  forging.  This  nonuniform  plastic  deformation  produces  residual  stresses  that  may  have  impor¬ 
tant  effects  in  the  strength  of  the  final  product.  Although  plastic  strains  are  unavoidable  in  large  deformation 
processes,  its  distribution  may  be  controlled  to  obtain  favorable  residual  stress  patterns,  in  particular  dur¬ 
ing  multi-step  forging.  Predicting  residual  stresses  calls  for  constitutive  models  that  account  for  the  elastic 
response  of  the  material.  For  this  purpose,  a  number  of  elasto-viscoplastic  models  that  use  hyper-  and  hypo- 
elasticity  together  with  internal  variable  representation  will  be  investigated.  In  particular,  constitutive  laws 
with  hardening  rules  that  reflect  microstructure  evolution  [4,  7,  24]  will  be  studied. 

Explicit  modeling  of  microstructure  features  as  represented  by  slip  processes  that  operate  during  plastic 
flow  as  well  as  the  reorientation  of  crystals  (texture)  comprising  the  metal  will  also  be  addressed.  These 
features  of  the  microstructure  are  introduced  in  numerical  simulations  of  metal  forming  by  polycrystal  plasticity 
models  [6,  19].  By  inclusion  of  these  models  in  the  firamework,  we  will  able  to  predict  texture  evolution  in 
the  forged  product,  and  hence,  the  optimization  process  can  focus  on  producing  an  optimized  geometry  of  the 
final  product  with  desirable  texture  and  mechanical  properties. 

Figure  9  summarizes  the  extended  framework  to  be  developed  in  Phase  II  for  the  optimization-based  design 
of  multi-step  forging  processes.  This  framework  will  be  complemented  with  a  methodology  to  be  developed 
for  deriving  the  corresponding  sensitivity  equations  which  are  needed  for  the  sensitivity  analysis  of  the  forging 
process.  The  specific  constitutive  models  to  be  used  in  the  framework  will  be  chosen  in  coordination  with  the 
Air  Force. 


7  Conclusions 

In  this  work  (Pheise  I),  we  have  developed  a  methodology  for  solving  optimization-based  design  problems  in 
linear  and  non-linear  mechanics,  along  with  the  numerical  tool  using  our  computational  substrate  PDESolve. 
The  solution  of  a  number  of  optimization  problems  using  linear  elasticity  and  small-  and  large-deformation 
plasticity  has  shown  the  validity  of  the  methodology  as  well  as  the  flexibility  of  PDESolve  for  implementing 
complex  tensorial  differential  equations.  In  particular,  the  solution  of  the  preform  shape  design  problem  for  a 
one-step  forging  process  have  shown  the  applicability  of  the  present  development  to  optimization-based  forging 
design,  a  technological  area  of  great  interest  to  the  Air  Force. 

In  Phase  II,  we  will  built  on  the  results  obtained  in  Phase  I  to  extend/improve  the  methodology  and 
solution  procedures  to  optimize  the  preform  and  die  geometries  in  multi-step  forging  processes.  In  this  case, 
the  computational  tool  will  be  based  on  the  finite  element  method,  a  version  of  PDESolve  which  is  currently 
being  developed. 
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Figure  25:  Schematic  roadmap  illustrating  components  of  constitutive  framework  for  modeling  forging  pro¬ 
cesses. 
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Appendix 

In  this  appendix,  we  derive  the  weak  form  of  the  governing  equations  for  both  the  direct  and  the  sensitivity 
problems  that  define  the  large  deformation  of  an  elasto-viscoplastic  material.  This  form  will  be  used  with  the 
finite  element  version  of  PDESolve  which  is  currently  being  developed. 

The  weak  form  of  the  governing  equations  for  the  direct  problem  is  obtained  by  taking  the  weighted  average 
residual  of  the  eqs.(10)-(12)  over  the  volume  Vn  (updated  Lagrangian  description)  to  give 

f  P:VudV=  f  pNh-udV-\-  f  t-u  dA  (95) 

Jvn  Jvn 

where  u  is  a  weighting  function.  To  obtain  the  weak  form  of  the  sensitivity  equations,  we  first  write  eq.(35) 
referred  to  Bo  (total  Lagrangian  description)  as 

Vo-<SP-Vo-(P-<5L2’)  +  P-(Vo-(5L2’)  +  ^pob  =  0  (96) 


We  note  here  that  the  derivation  of  this  equation  makes  use  of  the  reference  domain  Br  shown  in  fig.2. 
As  mentioned  before,  this  domain  is  typically  used  when  dealing  with  preform  shape  design  problems.  The 
weighted  averaged  residual  of  eq.(96)  can  be  written  as 

f  ((5P-P.(5Lj):VoudP=  /  6(>oh-udV+  f  [P.(V.<5Lj)]-udy  +  f  [{SP -P -SLJ) -no] -ndA  (97) 

J  Vq  J  Vo  d  Vb  d  Aq 


The  last  term  of  this  equation  is  equivalent  to 


—  •  ^Lo 


•  F^^:(n(g>  n)t  +  (<5Fu  •F^^:(I  -  n(g)  n))t]  •udA 


(98) 


which  is  derived  using  the  traction  boundary  condition  t  =  cr*n,  its  variation  (5t  =  (^o*— <T-^L^)*n+(n-JL^*n)t 
and  the  relations  (see  fig.  2)  =  JFx  •  F^^  Ft  =  F^-Fij,  JLq  =  JF/j  •  P  =  det{Fu)<^  * 

ndA  =  dei(Fu)F~^  •  ndAo-  Therefore,  the  weak  form  of  the  sensitivity  equations  can  be  written  as 


f  {6P-P-6L^):VudV  =  [  Spoh^udV -\-  [  [P.(V-JLj)]-uc/y+ 

Jvo  Jvo  Jvo 

[  [5t  -  F^  •  SLo  •  F-':(n  0  n)t  -h  (^F^  •  F-^:(I  --  n  0  n))t]  -udA 
Ja^ 


(99) 


If  the  reference  configuration  Br  coincides  with  the  undeformed  configuration  of  the  body  Bq^  then  ¥r  =  1 
and,  hence,  JLq  =  0.  With  this  value,  eq.(97)  reduces  to 


f  6P:VndV=  f  Spoh-udV f  [5t  +  (JF^  •  F-^:(I  -  n0  n))t] -udA 
Jvo  Jvo  Ja^ 


(100) 


Equation  (100)  is  mainly  used  when  dealing  with  die  design  problems  while  keeping  the  preform  shape  fixed. 

Discretization  of  the  weak  form  for  the  direct  problem,  eq.(95),  using  standard  finite  element  techniques 
leads  to  a  set  of  nonlinear  algebraic  equations  for  the  nodal  values  of  u  which  is  solved  using  the  Newton 
method.  On  the  other  hand,  discretization  of  the  weak  form  for  the  sensitivity  problem,  eq.(99)  or  eq.(lOO), 
produces  a  system  of  linear  algebraic  equations  in  the  nodal  values  of  5u  which  are  solved  at  each  time  step 
for  each  design  variable  once  the  solution  of  the  direct  problem  is  obtained.  At  the  end  of  the  specified  time 
interval,  the  computed  solution  for  the  field  variables  and  their  variations  are  input  to  the  optimizer  in  order 
to  obtain  the  new  values  of  the  design  parameters.  This  procedure  is  then  repeated  until  values  of  the  design 
parameters  that  minimize  the  objective  function  within  a  prescribed  tolerance  are  found. 


